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ABSTRACT 

Since its original proposal in 196I, the Particle-and-Force (PAF) 
method for numerical fliiid dynamic studies has been improved considerably 
and tested in a variety of new circumstances. This report, consisting of 
two parts, discusses at length a number of properties of the method from 
conceptual, practical, and experimental points of view. 

Part I is directly concerned with the methodology, giving in detail 
the procedure as it is now applied. It also includes the results of ex¬ 
perimental calculations, the conclusions, and a discussion of extensions 
now being developed. Part H delves more deeply into the meaning of the 
particle representation of fluid dynamics through a close examination of 
some pertinent idealized computer e:5)eriments. 
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PART I 


THE METHOD AND ITS APPLICATIONS 
introduction 

The Particle-and-Porce method was first described in a report [1-1] 
written in 1961 . The discussions given in that report are, to a great 
extent, still pertinent. Principal exceptions include the method of 
neighbor choice (page 12 ) and the manner in which the force function is 
derived (pages l6-19). Some of the more important parts of that paper 
are repeated in Part I here; others are omitted or merely summarized. 

In the Introduction to the above mentioned report, a discussion is 
given of the various techniques which have commonly been used for com¬ 
puter studies of time-dependent fluid dynamics in several space dimen¬ 
sions. The principal conclusion from that discussion is that no method 
exists which is universally applicable; each method has limitations. 

This means that the method must be chosen in each case to fit the require- 
ments at hand. There remain, however, many situations which still cannot 
be conputed by any one technique, or in many cases even by an appropriate 
combination. The situation is shown in general terms in Table 1-1. 

A second PAF-method report £1-2] riiows by statistical analysis how 
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t.e PA. interpartlcl. force can oe to correepona to the foro of the 

equation of state. These results have been extended and generalised in 

this report. 

In any practical PAF-method calculations, it is essential to use a 
large, high-speed electronic computer. The exas 5 les presented in this 
^ vere obtained with an IM-7C94 center, and confl^atlonal out- 
jut was processed through a Str»berg-C».lson S0-«>20 Microfilm Becorder. 

Chapter 1 presents a detailed discussion of all aspects of the PAP 
methodology fron. both conceptual and practical points of view. Some ac¬ 
tual tests of the computing technique are described in Chapter 2, while 
Chapter 5 contans some discussion of future plans for development, gen¬ 
eralisation, and extension of the PAP technique. Five Appendixes have 
been Included; they contan a r»ber of adationa topics suppl^enting 
the remarks in the body of the report. 
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CHATTEK 1 


DESCRIPTION OF THE METHOD 


A. DIFFERENTIAL THEORY 

The particles vbose dynamics we follow are to r^resent the elements 
of a fluid. Insofar as possible, we have used the classical particle- 
dynamics theory to edd in the formulation of the method. At some point, 
however, a divergence from the classical theory will be necessary so 
that the dissipative effects in a real fluid can be represented. Our 
pasrticles are not molecviles whose internal energy is carried by velocity 
fluctuations; indeed, we expect that the velocity of a particle is to re¬ 
present the mean velocity of the finite mass of fluid it r^resents. The 
macroscopic kinetic energy of the fluid is to be exactly the kinetic 
energy of all the particles, so the Internal energy must be represented 
by an additional variable. If this variable is esgaressed as a function 
of particle positions alone, then only adiabatic motion can be rQ>re- 
sented. Conqpression and subsequent expansion would then return the set 
of particles to exactly their initial configuration with no dissipation. 
Thus a specied prescription is needed to describe variations of particle 
internal energy. 
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we coneiaer the ayh«nloe of a eet of partlele. aeacrlhed hy the fol- 
lowing noinenclat\ire: 

i, j s indexes describing the particle number, 
s mass of particle #Jy 
r s space coordinate of particle #j, 
u s velocity of particle 

f S force exerted by particle #i onto particle #j, 

i.3 

^ S r /r (a tmit vector pointing from particle #L 
ij ij id ^ particle ^d), 

M s m u. s momentum of particle #dy 

V = 1 m n .u 3 kinetic energy of particle #d, 

^d “ 2 d d d 

j = internal energy of particle #d* 

nooenclature will te introaucea as requirea. 
we coemenee by assuBihg that the particles are govemea by the equa- 

tions of motion 


' L 


( 1 - 1 ) 


( 1 - 2 ) 


—r— = n.. 

dt 0 


The summation over i modified hy the presence of * does not include 
the term i = d, and is further restricted to include onDy certain 
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neighbors of j as discussed later. Stmmation without * includes all par- 
tides in the system • 

Next, we assume that the force function can be divided into two parts, 

(1-3) 


“ ^i/ij 


where 




The first term in the force is to be associated in form with the equa- 
tion of state of the fluid; the second term is introduced to achieve dissi- 
pation in the same manner as the "artificial viscosity" of Von Neumann and 
Richtmyer [1-9], or for the purpose of including real viscous effects. 
Correspond ence with Fluid Menban-fne 

Th. correspoi^ence vith fluid MohauloB ccoieB through bu eraaiuatlou 
of the conservation laws In fonns appropriate to the nature of the con- 
tlnum to he represented, to this discussion, the coordinate system Is 
Cartesian, cylindrical coordinates are considered later In this cheater. 

a. Conservation of mass. This Is automatic. Each particle has con¬ 
stant mass, mj, so that the total does not change with time. likewise 

the change of mass In any fired volume exactly equals the amount flowing 
over the bounding surface# 

Conservation of momentum. To satisfy this requirement, the re- 
rtriction is the same as in classical particle dynamics, namely 

% “ “^Ji* ^ demonstrate this, consider the momentum change rate of a 
particular subset of all the particles 
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Y m 
dt Zj J 

j(subset) 


I I V 

j(subset) i 


We break the sum over i into two parts and write 


i- I 1 , 

i(inside) subset) 


^ I , I 

l(outside) j(subset) 


where inside and outside refer to inclusion or exclusion from the subset. 

In the first double sum, each pair of particles enters twice, so that 

the total contribution to the sum from a particular pair is 

Since there must be no contribution to momentum change from particles 

the subset, the sum of the two terms must vanish. The second 

double sum does not thereby vanish, since each pair enters only once. 

Thus with P = -Pa 4 f the momentum change of any subset of particles 
i^ 

arises only through external forces, as required. 

The restriction also means that 


(l-i^) 


«i4 “ 


c. Conservation of energy . Here we must make a break from the 
usual procedure in particle mechanics. In devising a reasonable approach, 
we will be able to establish some of the crucial parts of the technique. 
The basis for the energy discussion is that the rate of change of energy 
of a particle should be given by the rate that the other particles do work 
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on it. This nork rate is in tur« given ty the product of force by velocity. 
(To be properly eyametric, the velocity through which the work flux is car¬ 
ried from one particle to another must be the mean value of the two velocl- 
ties.) Thus ve write 






( 1 - 5 ) 


It foUxws -ttmt the total energy of an isolated system is 


conserved. 


« at I ^ “ I * V - 


•> J 1 

Where ecjuality to zero follows from separate vanishing of the sum of con- 
trihutions from each pair. Likewise the energy of the particles in any 
subset changes only through work done on them by external particles. 

Now, we already know that 

dK du 

dT = wr 

Combination of this with Eq. (1.5) can be arranged to give 


dJ 4 -I ^ ^ 

dF" “ 2 X 


(1-6) 


Or 


(1-7) 


dJ. 

_i 

dt 


Z * dr ^ r—I* 

f —ii + 1 > ? •(? - 3*). 

^ ^ 2 L ' 1 y 


B, FINIIE TIME IKTERVAIS 

In practice, the numerical computations must proceed throu^ a 
sequence of finite time advancements whose steps are of duration 6t. This 
is accomplished throu^ a replacement of Eqs. (l-2), and (l-6) hy 


-♦ n+1 -♦ n 



5t 

-♦ n+1 





5t 

j 



■♦n+1 


&t 




(1-8) 


(1-9) 


(I-IO) 


where 




This shows how the variahles for cycle ^+1 are obtained algebraically 
from those of cycle ^b. The choice of time-centering of the equations is 


Ratified as follows: 


Eauation (l-8) — At the beginning of the calculations 

a cycle, the only information avrilahle fo^ha ^rce 
calculation is that which pertains to the beginning of cycle #a. 
The force is thus labeled with index n. 


Equation (l-9) 
there is some 


— After calculation of the new velocity by Eq. (l-8), 
arbitrariness as to what velocity should be used to 
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detemine the new coordinate. It vas shown on page 18 of Ref. 

newly calculated velocity is preferred to that 

^ of the cycle. The argu- 

ment is based upon requirements of computational stability. 

Eqmtion (I-IO) — The ri^t side contains the average of the 
old and new velocities, a combination which is introduced to 
ass^e rigor^s energy conservation in the time-difference fom 
oi the equations. (Mass and momentum have likewise been con- 
served; proof of this is the same as that for the differential 

equations. ) To demonstrate energy conservation we start from 
tne Identity 


1 

2 


Kr O'- - =’;) • 

Thus, from Eq. (l-8), the change in kinetic energy of a particle is 


i 


^i^tion of this with Eq. (1-10) gives for the change in total 
particle energy, E , “—o wuaj. 

u 


■r'^ E’ 


( 1 - 11 ) 


result shows that the energy transferred 

magniSde to but 

si^ from that transferred from #j to #i, thus prov¬ 
ing the contention of conservation. 


C. DISSIPATION 

There is at least one other form alternative to Eq. (l-IO) which 
could be considered for the internal energy calculation. The total energy 
difference could be calculated in a form analogous to Eq. (l-ll), but 
without ai^ time-centering of the ri^t side. The result would still be 
conservative as long as the proper reciprocal symmetries were preserved. 
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Prom the new total energy of the cell, the new kinetic energy could then 

he subtracted, giving the new internal energy. 

The reason for the specific choice of the form of Eq. (l-10) fol¬ 
lows from the requirement of monotonic dissipation. The use of Eq.. (I-IO) 
(with the proper choice of the g^^ forces) can never result (at least to 
lowest order in 6t) in decreased entropy, while most of the alternative 
forms can produce such decreases under certain circumstances. 

D. NEIGHBORS 

The manner in which a particle's nei^hors are chosen has evolved 
partly from the demands of the method and partly throu^ e^eriment. We 
already have mentioned some basic considerations of the method which place 
two primary restrictions on the choice of neighbors. 

1. Particles interact only with adjacent particles. 

2. A particle's effect on its neighlaor must be reciprocated, i.e., if i 
is a nei^bor of J, then J must also be one of i's nei^bors. 

The first condition is necessary for the representation of a fluid in which 

interactions are only with adjacent elanents, while the second is required 

for rigorous momentum and energy conservation. 

The nei^bor search technique which will be described here satisfies 
both of these criteria and, in addition, is particularly suited to effi¬ 
cient memory storage assignment. A search radius, R, and a maximum number 
of nei^bors, N*, are predetermined for each problem. Among all of the 
particles in the system which lie within the search radius of particle #j, 
the N* closest particles are tentatively chosen to be nei^bors. Any of 
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these ehlch do not slsmexa^r find J as one of their H* closest neighbors 
ere then subtracted from the list, thereby Insuring reciprocity. Good re. 

suits have been attained by choosing N* enual to twice the number of dl- 
mensions of the problem, 

?y restricting neighbors to only the / closest particles, this pro¬ 
cedure reduces the danger of interactions vith nonadjacent particles in 
highly compressed regions, where numerous particles Ue within the search 
radius. At the same time, however, the search radius restriction allows 
an interaction cutoff, which also can be useful. For instance, if one 
wishes to drive a shock throu^ a fluid at rest, R can be chosen less than 
the interparticle ^cing in the undisturbed material, so that this mater¬ 
ial will remain undisturbed until the shock reaches it. Also, when a par¬ 
ticle becomes detached from the main body of the fluid, it will cease to 
accelerate or to change its internal energy. This would not be true if 
there were no limit on R. In that case, the particle would continue to 
cool until its internal energy became negative and its force attractive. 

This meaningless local occurrence could play havoc with the rest of the 
calculation. 

An alternative method of achieving neighbor reciprocity would be to 
add rather than subtract neighbors, but the subtractive method has the 
distinct advantage for the programmer that he can place an upper limit 
on the storage requirements for interparticle relationships. Tests indi¬ 
cate that the two methods produce equally good results. 

The success of the PAP method d^ends to such a large extent upon 
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the pr<>per statistical averaging of interparticle fluctuations that It 
seems necessary to search for nelghhors every time cycle. It wouM he 
preferable to avoid this, because It Is by far the most time-consuming 
phase of the calculation, requiring 50 percent or more of the calcula¬ 
tion time. Therefore, Kg.eriment3 are being performed to see under 
.hat circumstances the neighbor search can be conducted less often. In 
addition, an attempt has been made to speed up this part of the calcula- 

tion as much as possible. 

1 ^ The Nel^hor Search Technique 

The first step In the current neither search technique is to over¬ 
lay the computational system vlth a grid of square cells, each ceU hav¬ 
ing a side of length R, the search radius. The particles are then clas¬ 
sified according to the cells in rtiich they fall and the actual search 
for neisibors begins. For any given particle i, the distances to all 
other particles vithin Its ovm or a neighboring cell are computed and 
compared .1th R. If any such distance, r^^, is less than R, then 1 is 
listed as a neighbor of 3 and J Is listed as a neighbor of 1, unless 
either 1 or 3 already has a full quota of neighbors. Then r^^ mould be 
compared to the distances to the other neighbors, and the most distant 

particle vould be axopped as a neighbor. 

m order to achieve the greatest possible speed vith this method, 
it is Imperative to choose R as small as accuracy viU permit. The cus 
tomary choice is about one and a half times the anticipated particle 
separation in the least compressed region. 
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If one starts at the loser left-hand comer of the cell mesh and 
works from left to right and upward through the cells. It Is possible to 
restrict somewhat the number of cells which must be searched for neigh¬ 
bors of any given particle. For example. If particle j lies in cell 

then It Is necessary to search only through cells (k, 1), (k,l + l), 
(k + 1, 1-1), (k + 1,1), and (k + 1,1 + i) for neighbors of if j 
had additional nel^bors In other cells bordering (k,l), this fact would 
have already been detemlned In searching for neighbors of those particles. 
After these tentative neighbor determinations have been made, a 

second pass Is made through the particles, dropping nel^bors as necessaiy 
for reciprocity, 

E. FCECES 

A theoiy Of correspondence between fluid dynamics and the pgp model 
has been worked out [1-2] which makes use of the methods of statistical 
mechanics. The conclusion Is that If the method can be developed In such 
a manner that each particle follows closely the Iccal mean of motion with¬ 
out extreme fluctuations, and If the force functions are chosen _ 

to a prescribed format, then the results should satlstically represent 
the fluid dynamics as desired, 

1• Equation of State Force 

The correspondence analysis [ 1 - 2 ] Indicates that the nomilsslpatlve 
part of the force function, fij(r„,Ji,Jj), should be associated with the 
equation of state of the fluid through an Integral equation which In two- 
dimensioiaal Cartesian ^ace has the form 
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( 1 - 12 ) 



•where 

m = particle density (assumed constant), 
p = particle density (particles per unit area), 
I = ^ecific internal energy, 
p = pressure (force per •unit distance). 


X = 

a as 


a. radial function describing the density of nei^bors, 
as defined in the previous section. The e^ct form o 
ij is not viniqueOy determined, the only restriction 
beiruz that it satisfy 


jt 



* 

d| = H . 


( 1 - 15 ) 


In Ref. 1-2, Eq. (l-12) is solved for the case in which o is approxi¬ 
mated by a step function. This formulation has proven satisfactory when 
applied to a po3ytrapic equation of state,* but for more complicated 
equations of state it may lead to a force function that does not vanish 
at normal density and zero internal energy. This difficulty can be 
avoided by choosing a form of a which is consistent with these require¬ 
ments. We now describe a method for determining such a form. 

Assume that the equation of state can be expanded in powers of the 


compression minus one: 


"^It can, in fact, be shown that every form of a which satisfies 
Sll^duce the same force function for a polytropic equation of state. 
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Ic ^ 

The nondlsslpatlve part of the force fpnctioa vlU then hare a ■h.-i.. 
appearance 

f(x,D = .“ ,J, 

k 

and the right-hand aide of E,. (1-1 a) yoXL conelet of a eerlee of inta- 
grals, one for each value of k* Since 


p(pq,I) = Aq(i), 

a must he of such a form that each such integral vanishes when 
except the one for vhich k = 0. 

Consider a typical term from the series. 


“P = Pq, 




I = px , to ohtfiiin 


<«/' 


Ct(-2 ^ m 2 

^ - 1 J CT(px^) dX. 


- ij »(e) at. 


®(l) d|. 


If the integrand is independent of n, then ^rCp^,!) » o, as it should. 
This suggests that 


p(l) * (constant) 8(| - i)^ 
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^er. 6 signifies th. Dir« f-ctlon. E^tlon (I-IJ) w .e u«d to 

evaluate the constant and ve obtain 
* 

o(|) *“ 7“ ^ 

With tbl. foon or 0, the correepondeooe e<suatioh hecooee 


p(oo,i) - 


Again, let I ® » 


- 5 d’^)- 


pUne coordloat... P - l/rjy .0 v. c«. finely «it. the cotreepon- 
aence heteeen the e.u.tloh-of-state pree^ sod the aohdiselpativ. pert 


of the force as 


f(rjyl) - sQ ’l) 


(1-14) 


Appljrlhg 


this fonmila to a polytropic equation of state. 


p = (7 - l)moI, gives 


f(rie,l) » -# (r - 1) -f- I 

H r. a 
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which, as indicated, is the same expression one would get using &ay form 
of o which satisfies Eq_, When we use Eq. (I- 14 ) for a "stif¬ 

fened" gas equation of state, p = (7 - l)tnpl + a(inp - p^), we obtain 



(1-15) 


Notice that the bracket term vanishes at normal density. Using for a the 
step function eB5)loyed in Kef. 1 - 2 , one arrives at the same expression as 
in Eq. (I-I5), excQ)t that the bracket term is 



which does not vanish at noimal density. 

The same corre^ndence result is derived by a virial approach in 
Appendix A, and it is also shown there why the number of neighbors per 
particle must be twice the number of space dimensions. 


Dissipative Force 


Reference 1-2 also provides a correlation between the dissipative 
force and true viscosity, so that theoretically one could innind** real 
viscous effects in the calculation. As a practical matter, however, the 
main requirement of a dissipative mechanism in PAP is, as in all finite 
difference methods which are not inherently stable, to prevent the growth 
of instabilities. In all of the e:!q>eriments reported here only the mini¬ 
mum amount of this "artificial" viscosity required for stability has been 
employed. 

*Ap|e^'I.°^ correspondence for dissipative forces is discussed in 
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Dissipation is only required in those portions of the fluid which 
are undergoing coinparison; hence the dissipative force, g^j, between 
particles i and J is nonzero only when - u^) > 0. Ih a previous 

report on the PAP technique [I-l], a dissipative force of the type 

hi ■ "“<^1 - 

was proposed, where co was a constant with dimensions of reciprocal time. 
But, in Appendix II of Ref. 1-1, it was demonstrated that this form of 
the dissipative force did not rigorously conserve angular momentum; there¬ 
fore, an alternative of the form 

-» r“4 /-^ i? 11 (I“178-) 

ms suggested. Hot only does a force of the type Eg,. (l-17a) conserve 
-C-t-r »n.entu», hut it also gives a measure of the compression hetwen 
particles i and J. [Hote that Eg.(l-17a) represents the coaiponeiit of 
Eg. (I-l6) lylee along the Une of particle centers.] Thus Eg. <I-17a) 
is consistent vith the oncdlmensional stability analysis uhlch assumes 
that all compression is head-oni it is Uss Utelly to lead to interpene¬ 
tration than is Eg. (1-1 6 ) for a large (but theoretically stable) value 

of a>. 

The two forms, however, share the drawback of being independent of 
the interparticle ^cing. Thus the dissipative force is the same when 
the particles are a search radius apart as when they almost coincide. 
This situation has been corrected in subsequent work by ejaploylng one or 

modifications of the type Eq. (1-17a) force. 


the other of the following 
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where <a is now diaienslonless and is chosen equal to a typical sound 
speed for the prohlem. The first of these two forms has been in the most 
common usage; the second one is reserved mainOy for problems whose initial 
conditions require the material to be cold. The force Eq. (l-17b) is 
patterned after a form suggested by landshoff [ 1 - 10 ]; the square-root fac¬ 
tor is simply proportional to local sound Epeed in a poiytrople gas and 
closely related to it for many other materials. 

P. CYUNDRICiLL CQCtRDTKA»PRa 

Many fluid dynamics problems deal with flows which possess an awgni a-i* 
sysanetry about some axis throu^ the flow field. To allow the PAP method 
to take advantage of this j^mmetry requires certain modifications of its 
concepts. The particles now become rings about the axis of (^vietxy and^ 
since rings at different radii occupy different volumes, it has been con¬ 
venient to make the mass of a particle proportional to its -tMi t lal distance 
from the axis if the d^sity throughout the fluid is constant. 

At the present time, no correlation theory such as that given in 
Ref. 1-2 exists for cylindrical coordinates. (Some newly derived results 
are given in /q>pendlx C.) Therefoire, in order to arrive at a force func¬ 
tion appropriate to these rings about the axis, it was necessaxy to resort 
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to the equations of motion themselves* This approach lacks the unam- 
hlgulty of the theoretical correlation method, so there were, as might he 
eacpected, several false starts. We review the major one of these in 
order to demonstrate some of the rather subtle difficulties which can 
arise with PAP and to show how these were eventualiy avoided. 

1 • A ft^HTnitiayy CyH^^'>*ical Coordinate Force 

In flows possessing cylindrical symmetry, a c ylindr ical shell of 
fluid is seen to have a radial acceleration which is independent of 
exterior forces and d^nds only r^n the pressures that exist within 
the 1 (for example, see pages 55 3^ Ref. 1-11). The average 

radle^. acceleration of such a ^ell is given hy 
du 

® iSr* * ^^^^inslde ” outside 

where A refers to the surface area (inside 
thickness and length L. If we 

consider the PAF particles to he arranged 
in a rectangular array (as in the acccn- 
pemylng figure),then, since each particle 
actually represents a toroidal shell of 
fluid about the cylindrical axis, Eq. 

(l-l8) might he suitable for the derivation of the radial nondissipative 
forces felt hy a particle. Assuming a polytropic equation of state, 
p a (7 - l)pl; ( 1 - 18 ) can he rewritten 


2)tL 


/: 


P ar. 


( 1 - 18 ) 


or outside) of the shell of 
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vbere v (= u^) is the radial component of velocity and L is taken to be 
the inteiparticle gjyaclng in the z direction. Since the treatment in 
the z direction is entirely unchanged from that of plane coordinates, it 

would seem that the nondissipative portion of the force felt by particle 
j could be written 



where r is a unit radial vector. 

The self-administered force, the second term on the right in 
Eq- (1-19), is the force exerted on the particle by the rest of its own 
ring. It is positive for positive J and can be thou^t of as a self- 

ejq>andlng force, tending to carry the particle away from the axis of 
symmetry. 

Extending this same meaning to the dissipative part of the force, 
one can view the self-administered part as the conpressive force which re- 
suits whenever the radius of the ring decreases. The force is in the 
positive r direction, so we could think of it as the force exerted by 
the image of the particle through the axis of symmetry. Therefore, the 
self-administered part of the dissipative force should be obtainable from 
the normal dissipative force with the velocity of the nei^bor, re- 
placed by -u^. At the time that these forces were derived for cylindrical 
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coordinates, a dissipative force vhlcli was the sum of the e::^ressions in 
Eqs. (1-1 6 ) and (1-17a) was heing used. Hence, the self-administered dis¬ 
sipative force was taken as 

^ j~2a^(ai^ + <*> 2 )Vjr when 5 ® 

^0 when v > 0, 

J 

Rnfl the entire force exerted on particle j was 


=r{(^ 


(j, * j.) 








- 1)J, 




( 1 - 20 ) 


for a polytropic equation of state. Here, m^^ = l/2 (m^ + m^). 

Several problems were run using this force function. When applied 
to a problem involving flow past a cone, it was dlscoveired that the shock 
reflected from the face of the cone did not take the form expected, but 
became elongated with time in the radial direction. This difficulty was 
traced to the fact that the self-ejqoanding force was producing an exces¬ 
sive mass flux away frcmi the axis of symmetry. 

As a result of this dlscoveiy, the derivation of Eq. (l-20) was re¬ 
examined. It was reallz€$d then that a cylindrical shell of fluid bounded 
by flree surfaces would expand both Inward and outward at escape velocity, 
and that the net increase in outward momentum was simply a result of the 
fact that more mass was e:q>osed to the outer surface. The nonnegative 
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self-e 3 cpandlng force was therefore clearly wrong. What was needed was a 
force function which wotild eUJow an unconflned fluid to es^and In an 
directions in a manner that was consistent with the equations of motion 
in cylindrical coo3rdlnates. 

Very often in the derivation of finite difference formulas it is 
found that equivalent differential es^resslons can give rise to dlffe 3 rence 
equations with markedly different properties. Such is the case here; for 
if we return to the momentum equation with cylindrical symmetry, 


dif 9p 
Pdt“ 



( 1 - 21 ) 


from which £q. (l-1d) was ohtedLned, we can derive force esipressions which 
do not give rise to these large mass fluxes. 

2. The Rresent Cylindrical Coordimte Force Function 

Once again, let us consider the particles to he in rectangular array 
and examine the effect of the r component of Eq. (1-21) along a radial 
line in this array. By referring to the particles along this line by an 
index number, k (which increases with radius), we can write 


1 

^ TIT - to M ■ W’ 


where we have added curblflclal viscosity tezms, q, in the form of pres¬ 
sure modifications. Assume that 
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if > 0 
otherwise;) 


r ^‘“ Pk ± iW’'w ■’' k ’ 

s±i° to 


where the same as yet luidesignated quantity having the dimension 

of a velocity. Then we can rewrite Eq. (l-22) as 




(r-i)r *^k 

'^k-l 

"^k+l 


_1 



dt 
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in Eq. (I-I7b), ^ 

we choose 
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qrtfi maXe the substitution 
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2r. 


^k+1 " 

this equation can be expressed as 


"Vat 


(7 - Dxj. 


r Jk * Jk-1 _ 

• ’'k-1 ' ^+1 






r^vi - \>y 

+ 

1 

+ ojrj^ 





T 


" \+i^ J~ 

"'k^t ^ ^k 

^ VI 

(vi * ’^k' 




-55- 




Thus, generalizing to any random array of particles, we obtain for the 
r component of the interparticle force the expression 






(l-23a) 


^ere the term with m in it is set to zero if not positive. The 
2 component of the force carries over directly from plane coordinates. 



Notice that this choice of produced a g^^ force of the type Eq. 

Cl“17b). If we had chosen s constant, we would have obtained a force 
of the type Eq. (l-17c). 

A study of the function e^qjressed in Eqs. (I- 25 ) confirms that it 
does satisfy the requirements demanded of a force function for calcula¬ 
tion in cylindrical coordinates. First of all, having been derived from 

the equations of motion in that coordinate system, it is consistent with 
them to at least first order. 

The other quality which we demanded of a cylindrical coordinate 
force function was that it allow an unconflned fluid of uniform density 
and internal energy to have the same rate of ej^wmsion in the positive 
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and negative radial directions. The force function e:<pressed in E^s. 
(j_ 25 ) does have this property. To see this, consider a system in vhi h 
all particles have specific internal energy, Iq. The radial acceleration 
in^jarted to particle hy its neighbor, #i, is 



in this uniform density fluid, particles are equally spaced and have 
masses proportional to their radii. All particles, therefore, e^erience 
the sane acceleration, so rarefactions move into the fluid from the inte¬ 
rior and exterior boundaries at the sane velocity. 

conservation of energy in a cylindrical coordinate system requires, 

as it did in the plane coordinate case, that the rate of change of the 
energy of particle j be given by the rate that other particles do vork on 
it. This vork rate is in turn given by the product of interparticle 
forces tines appropriate velocities less the portion of this product which 
is attributable to self-expanding forces. Since the total energy of the 
system nust be conserved, we have 

^ ^ ^i ^ ^ (Energy of Self-Expansion)^ = 0. 

Thus the rate of change of the energy of self-expansion is given by 



where [ 1 represents the bracket term in Eqs. (l-25)j 

4.J 
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so we can write the eqxxation for the rate of change of internal energy in 
cylindrical coordinates as 


dJ dE. dK. 

- ^ =: _J - _ ± 

dt dt dt 


‘ i M’j'V'ij ’ 


i I. -'ij ^ ^ 


■“ji 


( 1 - 24 ) 


In the finite time-difference approximation to Eq. (1-24), the velocities must 
he time-centered for energy conservation, as was the case in plane coordinates. 

By extending the procedure used to derive the polytropic equation of 
state force, Eqs. (1-23), to the "stiffened" gas equation of state, one can 
derive the analogue of Eq. ( 1 - 15 ) m cylindrical coordinates. 






r 




.2 1 


(7 - 


1 ) (Jj + J.) + a(in + m ) fl-- .l^^ij^o l 


(l- 25 a) 
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+* •'ill"! * ”j) j 

(7 - 1)(Ji + J ) + a(m + m.) fl - — 1 
^ J i J L (m^ + mj) J 


(l-25b) 
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This force function has been applied successfully to several picblems, one 
of vhich is described in Chapter 2. 

G. BOOMDARY COKDITIOMS 

From a computationeQ. standpoint it vould be perfectly reasonable to 
perform a PAP calculation vithout boundary constraints, since the method 
is not tied to a fixed mesh system4 But for the practical reasons of effi¬ 
cient neither searching, taMng advantage of flow symmetries and the log¬ 
ical presentation of results, the fluid is usually bounded. In many phys- 
circumstances, of course, such bounding also exists. 

Two types of boundary conditions have been used (sometimes in com¬ 
bination): prescribed input, continuative output conditions for flow into 
and out of a channel, and reflection from rigid walls. 

1 , Prescribed Input, Continuative Output 

In the case of flow through a channel, particles are periodically 
cjjeated at the input boundasjy with the proper mass, velocity, internal 
energy, and interparticle spacing to represent the required state of the 
incoming fluid. At the output boundary, particles are destroyed and 
their mass energy are subtracted from the totals for the system. The 
memory storage which had been used by the destroyed particle is then made 
available for new paxticles, so the problem can be run indefinitely. 

Particle creation and destruction actually occurs at a distance 
of approximately the search radius outside the input and output boundaries, 
zmspectively. During the time that these petitlcles exist outside the 
system, they exert a force on their neighbors inside, but they themselves 
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are not so influenced. They therefore do work on the interior particles, 
so it is necessary (for the purpose of checking energy conservation) to 
calculate the change in the energy of the system vhlch stems from their 
influence, m plane coordinates, the change in energy ^ch results from 

particle outside the system exerting a force on its nel^hors i in side 
the system is 


5£ s 5t 




. 'i'- 


(1-26) 


But in cylindrical coordinates a portion of the inteiparticle reaction is 
attributable to the self-eaq>anding force of a particle. This is not 
actually vork done by an external body on the system, so the energy change 
in Eq. (I- 26 ) must be reduced by this amount. The self-expanding energy 

attributable to each of the neighbors, i, above, is given by one-half the 
amount by which 


r - ^ ^ (i?... ?,)■, 


fail, to reduce to eero. Heoce, l,y maJilng uee of the fact that in cyllu. 
drlcal coordinates 

we can write the cylindrical coordinate analogue of Eq. (I- 26 ) as 


SE ss 5t 


r { fe- - i c^) c-^)}. 


(1-27) 


The velocities in Eqs. (l-26) and(1-27) must he time-centered in the 
finite time-difference form, as was required in the previous sections on 
energy consezvation* 

2* Rigid Walls 

The problem of devising boundary conditions which accurately describe 
the interaction of a fluid with a rigid wall is one which has received 
extensive study. We have experimented with many different boundary treat¬ 
ments, almost all of which work qpite well in plane coordinates or in un¬ 
obstructed flow parallel to the axis of syrametiy in cylindrical coor¬ 
dinates. The most difficult problem for all such methods involves flow 
past a cylindrically symmetric obstacle in a channel. 

The difficulty arises from the fact that the boimdary conditions at 
the obstacle cause rather large fluctxiations in what shoxild be a stagnant 
region adjacent to it. These fluctuations are compounded in cylindrical 
coordinates by the presence of both very light and heavy particles at the 
face of the obstacle. A light particle caught between the wall and several 
heavy neighbors may begin to oscillate at such a rapid rate that it becomes 
drained of internal energy and develops a negative (attractive) force. In 
time, this particle’s neighbors will also develop negative forces and the 
results will soon become meaniisgless. 

We have atteurpted to alleviate this tendency by developing boundary 
conditions which damp particle fluctuations near the wall, but in some 
cases more drastic aotion is required. Fcr example, c6p.culations of flow 
past a cylindrically symmetric cone weie attempted using each of the 
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boundary conditions described In the next two sections of this report, and 
In no case couOd the problem be run long enough to achieve steady state. 

As light particles moved up the face of the cone and encountered massive 
neighbors, boundary Instabilities Invariably developed. One way of avoid¬ 
ing this difficulty Is to destroy the Ugjht particle before the instability 
fluctuations damage the calculations# 

Any particle whose mass is less than some fraction, k, of the average 
mass of its neighbors Is destroyed, and its mass, moment\aii, and energy 
are distributed among its neighbors. Values of k from 0.2 to 0.4 have 
been suitable for this purpose* 

The mass of the destroyed particle is distributed among Its neighbors 
on the basis of interparticle distance. That is, a particular neighbor, fL, 
of the light particle, #J, gets some fraction. 



of J*s mass. The new mass of particle §L Is 


“i = “l + 

where the tilde indicates the quantity before distribution. 

The total momentum and energy of the group of particles will -fe h o n be 
conserved If we canpute each neighbor’s new velocity and Internal energy 
from the formulas 
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tojther proceaure, descriiea in y^peniin E, nnkes use of a mass 
flux from heavy to light particles. 

5 , Mecbaplcal Reflection fr om a Wall 

Ivo general classifications of DounOaxy conditiOTis, along vitn manir 
vsnations. hare teen enployed: nechanical reflection of particles fro. 
tbe wall and reflection ty neans of image particles on the opposite side 
of the wall. The simplest type of mechanical reflection (from a pro- 
gra.«r's standpoint) is to allow the particle to cross the wall, calcu¬ 
late its interpartlcl. forces, its weloclty, and its IntemaL energy in 
the usual imnner, and then, for the purpose of moving the particle, re¬ 
verse its component of velocity normal to the wall. Energy is conserved 
hy this process and, sine, tha particle will he accelerating as it 
crosses the houndary, it will re-enter the system within one time cycle. 
But, rather than having a dissipative effect on the fluctuations of 
adjacent particles, this manner of treatment increases those flnctua- 
tlais hy returning the particle to the system with a higher velocity 
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than that with which It left. To alleviate this sitmtion, an attearpt 
was made to allow the wall to exert a force on nearby particles within 
the system, thereby slowing them down before they crossed the boundary. 
This, however, had the disadvantage of subjecting particles near the 
boundary to continual compression, thereby raising their internal energy 
to abnoamially high levels. 

Aa alternative approach toward reducing fluctuations near an obstacle 
wall is to set to zero the normal component of the velocity of any par¬ 
ticle which strikes the wall. The particle is therefore constrained to 
move in the plane of the wall until it reaches the top of the obstacle. 

The Physical basis for this boundary condition is the fact that the nor¬ 
mal component of velocity of a continuous fluid impinging on a rigid 
wall would be zero. 

But, while this approach danps out flact\iations normal to the wall, 
it leads to severe tangential fluctuations. Since the particles at the 
wall are very hot (the lost kinetic energy is transformed to internal 
energy) and highly concentrated (since none can escape from the wall), 
their interparticle forces become extremely large. The great accelera¬ 
tions which result will cause particle Interpenetrations unless the time 
step is made very small. 

The most successful form of mechanical reflection of particles from 
a wall has been specular reflection. At the Instant that a particle 
strikes a wall its velocity component normal to the wall is reversed. 
Therefore, there is no time lag in returning a particle to the system. 
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as there vas aith the treatment hy (*leh a parttele actually penetrated 
the uall before being reflected. This greatly reduces the fluctuations 
in the vicinity of the vail. Icicles tend to sticl to the sail uith 
this method, as they did in the case idiere the normal velocity component 
vas set to aero. Mow, however, this is not so serious a problem, since 
the particles on the wall ere neither especially hot nor highly concen- 
trated, as they vere in that case# 

4. Bnages 

All of the forms of mechanical reflection from a vail share one im¬ 
portant defect: they affect particles in a vay that is entirely unre¬ 
lated to the interparticle relationships which exist in other parts of 
the computational system. Granted that rigid walls must, hy definition, 
present an abrupt obstruction to the flow, it is usual3y preferable in a 
computational scheme to make this treatment an extension of the calcula- 
tional method in the interior of the system. This can be accomplished 
vlth PAF by the use of self-images for particles near the boundary. 

An image particle is created (see the following diagram) for every 
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particle In the cystem which lies within a distance of one-half the search 
radius of a rigid boundary and from which a normal to that boundary can 
be constructed. Notice that with this definition, it is possible for a 
particle within the required distance of a wall to have no, one, or two 
images* A particle at a concave boundary corner would require two 
images so that it wouM be aware of both boundaries. At a convex boun¬ 
dary comer it would be possible to create an image which is a reflection 
of the particle throu^ the point of intersection of the two boundary 
lines. This was not done for the problems described in this report, but 
subsequent experimentation indicates that it is desirable. 

The image particle is given the reflected properties of the corres¬ 
ponding particle within the system. If particle i is the image of a par¬ 
ticle i within the system, these properties can be described as follows. 
Coordinates: Particle i is given the coordinates of a point ^ch is the 
same distance from the wall as j and lies on the same normal to the wall. 
These coordinates can be written: 



f 

where 5 is the slope of the wall and x^ is its intercut with the x axis. 
Velocities: Particle i has the same component of velocity tangential to 
the wall as does J, but it has a normal component which is the negative 




Mass: ®i ** 

Internal energy: 

Of course, these properties apply only to straigtht boundaries. To devise 
boundary conditions that are appropriate to curved boundaries vill require 

additional thought and e^erimentation. 

B>ages are created only when needed and are destroyed after each cycle 
of calculation. The fact that a particle is sufficiently close to a rigid 
wall to have an image does not necessarily mean that the image will be one 
of the particle's nei^bors. If the particle should find N* closer neigh¬ 
bors, then its image will have no effect. On the other hand, it voiild be 
possible to include among a particle's neighbors the Images of other par¬ 
ticles. In this case we would have to require (for energy conservation) 
that an additional principle of reciprocity be adhered to, i.e., if i 
finds ya image as a nei^bor, then J must have i's image as one of its 
neighbors. None of the problems described in this report involve such 
multiple image neighbors. 


5 , Conservation of Energy 

linages have as their sole nuipose for existence the influence of 



other particles; hence, we neither calculate changes in their properties 
not include their mass and energy in the totals for the system. Never¬ 
theless, they behave in every way as conventional particles and, for the 
purpose of confirming the presence of energy conservation near a wall, ii 
is convenient to consider them as such. 

To demonstrate that the image method 
req,ui3res no special treatment to insure 
consesrvation, let us consider a specific 
exas^le. The accompanyiiag figure shows 
two particles, 1 and 2, and their re¬ 
spective images, 1 • and 2', across a wall 
in either plane or cylindrical coordinates. 

no effect on any of the calculations, we are free to look on these four 
particles as an isolated region of the fluid and allow them to Interact 
only with one another. Total energy is conserved among these four 
particles, as it would be among any isolated subset of particles, i.e., 

6E^ + SEg + 8E^, + SEg, = 0. 

Furthermore, by symmetry, the work done on any one of the particles 
during a time cycle (vhich equals the change in its energy that 

cycle) is equivalent to that experienced by its image; hence 

SE^ + SEg = SE^, + bEg,. 

To be consistent, the above equations require that 



• I 

I 


Ince the wall itself has 
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5E^ + BEg ss 0, 

i,e*, that the energy of the interior particles he conseinred. 

Because of the asymmetry of the r component of the force function 
in cylindrical coordinates, it mi^ seem that the work done on the 
particle and its Image would not be the same in that coordinate system. 
But, in cylindrical coordinates, the actual work which a particle i does 
on its neighbor J is defined (for the PAP method) to be the plane coor¬ 
dinate work rate. 



less the pcwtion 




* V 



which is attributable to the work the ring i does on itself. Thus the 
actual work rate is 


^(F 
2 ^ i^ 



and this is a symmetric function. 

H. PARTICIE COMFIGURATIONS 

A problem was run rather early in the development of the two- 
dimensional PAP technique in order to determine the equilibrium con¬ 
figuration of a box of particles. The box was initially filled with 
a random arrangement of warm particles at rest. These particles then 
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began to interact and it was observed that \dien the initial large fluctua¬ 
tions had subsided the particles had attained a hexagonal configuration. 

This sane hexagonal array of 
particles (see acconipanylng figure) 


was observed in every problem in 
which a rectangular (or near rectan¬ 
gular) configuration was subjected 
to some type of defoming process. 
For example, as a shock moves 
through a rectangulsir array of 
particles at rest, such a configura¬ 



tion is attained within two or three particle i^acings behind the shock. 
If the particles at rest are initially slightly displaced in some random 


manner, the hexagonal array is assumed ImmediateOy behind the shock. In 
this case the shock is observed to move at more nearly the theoretically 

predicted speed than in the case where the particles are not subjected to 
a random displacement* 


The primary advantage of random particle displacement is that it 
allows systems of particles to become deformed rather easily. When a 
disturbance moves through a system of rectangularly arranged particles 
in a direction parallel to their rows, the two-dimensional effects of the 
distxurbance are resisted because the particles can exert forces on one 
another only along a line connecting their centers. Therefore, in almost 
all two-dimensional problems it is preferable to stagger particles 
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subtly from a purely rectangular array. The only notewortly exception 
to this rule occurs when, for illustrative purposes, one wishes to see 
where the first slight two-dimensional effects become apparent (as in de¬ 
fining the bow wave in flow past an obstacle). 

Particles are generally offset as much as one-fouarth the inter- 
particle spacing in each of the coordinate directions. Two random numbers 
between -1 and +1 are generated for each particle and the particle is then 
displaced from its rectangular position by these fractions of the maximum 
displacement in each of the coordinate directions. 
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CHAPCER 2 

TESTS OF THE METHOD 


A previous report on the PAP method [1-1J Included de«a^ptlons of 
several tests of the method. These vere simple fluid floe prohlems for 
ehleh anal^rtlc solutions existed. The calculations ears In good sgremsent 
With these analytic solutions. 

we nov describe some more complicated tests of the 
technique, for *lch analytic solutions are not available. The results 
Of most Of these are checked against «q«rlmental results from shock tube 
end wind tunnel tests, and one Is compared with the results obtained by 
another numerical method. Several other problems, for which there are no 
comparisons, are used to Illustrate some of the properties and capablU- 
ties of PdP calculations. The ccrngnitatlons vere performed on IBl 
Electronle Data ITocesslng Machines types 7090 and 7094. 

A. FIW PAST A WEDGE 

This is a study of the rate of growth of a detached how wave pro¬ 
duced by the passage of a shock over a two-dimensional wedge su^nded 
in air. Time zero corresponds to the time of arrival of the shock at the 
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apex of the wedge! we therefore heglh the PAP problem with parttclee In 
a shocked state on the left of a wertieal line through that point and In 
an ambient state on the right of that line. Particles move from left to 
right pest the wedge and exit through the contlnuatlye right-hand boundary. 
Sew particles. In a shocked state, enter the system throujdi the left-hand 

IxiuMazy* 

per this problem the Mach number of me flow is 1.55 and the «>ex 
angle of the wedge Is 90°. We use a polytropic equation of state for the 
air, with specific heat ratio r - 1.4. Other input data are listed below, 
where subscrlirt 1 signifies shocked air end subscript o Is ambient air. 


X component of velocity (cm/|isec) 

u, = 0.0755 

u = 0 

0 

y component of velocity (cm/nsec) 

Vi-0 

Vq = 0 

particle internal energy (10 joule) 

» 1.0403 

= 0.5952 

particle spacing in the x direction (' 

cm) Sx^ = 0.0655 

6x = 0.1245 
0 

particle spacing in the y direction ( 

cm) &y., = 0.0655 

6y^ = 0.1243 


/ N m = U86l X 10 

particle mass vgm; 


time st^ (nsec) 

artificial viscosity coefficient 

(type I-17a force) “ * 

The locations of the detached bow wave at times of 18 and 46 nsec 

after Impact are compared In Fig. 1-1 with those observed by Griffith 

tI-12) in a shock tube experiment. The PAP results are Indicated by the 

dashed lines In the figure, while the solid lines show the location of 

Griffith’s detached ^ock fronts at the tvo times. 



The shapes of the two curves are someiAat different at the first 
Observation tl«e, but are essentially the sene at the second tine, except 
for the reflection on the FAP curve ^ch was caused by Interaction with 
the top reflective boundaiy. Griffith reports, however, that his Measured 
ehoch strength of , .35 nay he In error by as nuch as %, so the guallta- 
tlve agreenent between the two curves Is as nuch as could he e^wected. 

This in Itself Is Inportant, for It demonstrates (at least qualitatively) 
that PAP calculates early bow wave development correctly. 



fig- I-l. A conparlson of the PAP detached bow wave 
positions (dashed lines) at 18 and 46 nsec 
after impact with those ohserved in a shock 
tu^ e:q)eriinent involving a Mach 1 .35 flow 
past a wedge. 
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The Griffith report shows subsequent bow front positions at times up 
to and including steady state, vhich is attained at about 5OO nsec, 
tfofortunately, computer memory limitations restrict the PAF run to about 

90 psec. 


B. FLOW PAST A BIIJMT CYLINDER 

In contrast to the previous problem, in which we were concerned 
with the ability of PAF to resolve the initial characteristics of flow 
past an obstacle, we now consider a test of its steady-state predictions. 

This is the study of a Mach 1 .58 air flow (7 = 1 .^) Past a blunt 
axially symmetric cylinder in cylindrical coordinates. We are not con¬ 
cerned with the time dependency of the flow in this case; hence we can 
state the free-stream conditions without specifying their units. Those 

conditions are: 


z coaponent of velocity 
r component of velocity 
sound speed 

particle spacing in the z direction 

particle spacing in the r direction 
psorticle mass 
time step 

sirtificial viscosity coefficient 
(type 1-17b force) 


u = 1.58 
o 

Vo=0 

Co= 1.0 
5Zq =1.0 

6 r = 1.0 
o 

m = initial radius 
8 t = 0.1 

<D = 0.5 


The first test of the calcxalations is a comparison of the PAF steady- 
state stagnation pressure with the theoretical value. Tne theoretical 
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value can be determined by making use of BesmouUi’s law, which states 
that along a streamline in steady flow, the quantity 



2 . 

+ V ) + 



P 

P 


is constant. We can therefore relate conditions at the stagnation point 
(u 3 0, V s o), which is located at the center of the cylinder face, to 
conditions at a point immediately behind the detached shock on the axis 
of symmetry, because these two points lie on the streamline. Referring 
to values at the stagnation point with subscript s and to the ahncked 
quantities with a prime, we can write this relationship as 



+ v'^) + 



The stagnation density can be eliminated by means of the adiabatic equa¬ 
tion of state. 



to obtain an expression for the stagnation pressure in terms of the con¬ 
ditions behind the detached shock. 




7/7-1 


(1-28) 


The shock conditions can, in turn, be i^lated to our input quan¬ 
tities by means of the Rankine-Hugoniot relations. The theoretical value 


of the stagnation pressure, Eq. ( 1 - 28 ), for the present problem is 

p = 0.426. In Pig. 1-2 the PAP pressures are compared with this value 
^s 

for times 15 through 37 « 

Although the PAP pressures fluctuate a great deal, the average of 
these fluctuations is in good agreement with the theoretical stagnation 
pressure. 

The reason for these fluctuations is that the pressure at any 
instant of time is obtained by summing the forces exerted by image 
particles across a small disk at the center of the cylinder at that time 
and dividing this sum by the area of the disk. The number of particle 
interactions sampled at a particular time is Small, so the total force 
varies considerably wltb ‘time* 

Another check on the PAP steady-state predictions can be obtained 
by comparing the final PAP bow wave position with that obtained by 
Marschner [1-13] for this same problem in a wind tunnel e:5)erlment. 

Figure 1-3 shows a steady-state PAP particle plot on which has been super¬ 
imposed Marschner's steady-state detached shock position. Good agreement 
is evident near the axis of symmetry where the shock is strong. In order 
to provide a basis for comparison in the weak shock region away from the 
axis, the PAP particles were constrained to enter the system through the 
left boundary in a linear fashion (in contrast to their usual staggered 
arrangement). The point at which these lines of particles become de¬ 
flected from their horizontal trajectories is the point of intersection 
with the PAP shock front. One may verify that this shock front agrees 
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stagnation pressure observed in a PAF calculation of a 


13 

ll 

o S 






Pig. I- 



The steady-state detached shock fJront position ohse^ed 

eiqpertmentaUy in a Mach 1 .g flow SSSe 
sywetric cylinder is plotted on a late-time PAF^mcie 
configuration from the calcvilation of the same problem. 
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nicely with Marschner’s bow wave. Because of the weakness of the shock 

in the vicinity of the reflective boundaries, there is no visible evidence 

of shock reflection as there was in the previous problem involving flow 
past a wedge. 


C. FLOW PAST A COKE 

The third and last example involving flow past an obstacle is con¬ 
cerned with the passage of a Mach 1.41 air shock over a 75° cone. We com. 
pare the steady-state pressure along the obstacle face and the final bow 
wave shape with those observed experimentally by Marschner [ 1 - 15 ]. 

The free stream conditions for this problem are; 


z component of velocity 
r component of velocity 
sound speed 

particle spacing in the z direction 
particle spacing in the r direction 
particle mass 
time step 

artificial viscosity coefficient 
(type 1 - 17 b force) 


u^=1.4l 

c =1,0 
o 

8z = 0.1 
o 

6r^ - 0.1 

m = initial radius 
5 t = 0.01 

<15 = 1.0 

the PAP pressxires along the face 


Figure 1-4 shows a con^iarison of 

of the cone (the dots) with those observed by Marschner (the line). 

The location on the cone face is plotted in the nondimenslonal fonn x/j 
where x is the distance along the face of the cone measured from the nc 
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and s Is the sleuifc length of the 


cone from nose to shoulder. The pres- 
sure, arc also plotted In n nondlmenslonal form p/p^, phere p Is the oh- 
served pressnre at the point and p^ Is the theoretical stagnation pressure 
vhlch develops on a blnnt cyUnder subject to these free stream conditions. 
The stagnation pressure used by terschner Is an observed value, vhereas 
that used to reduce the EAF data uas calculated from Eq. (i-aS). 

The ESP pressures were obtained In a manner similar to those plotted 

in Mg. 1-2, The z components of the forces enerted by Image particles 

across a punctured dish of radii x/s ± 0.1 were summed, averaged over time, 

and divided by the area of the punctured disk (this Is equivalent to 

dividing the total force across a segment of the cone face by the area of 
the segment). 

Agreement with Marschner's curve Is considered to be good, although 
there Is less variation In PAP pressures across the cone face than was 
Observed In the wind tunnel e:q«rlments. Perhaps this can be attributed 
to the PAP strategy (described In the discussion of boundary conditions) 

Of destroying light particles as they move .q, the cone face and distribu¬ 
ting their mass, momentum, and energy among their heavier neighbors. As 
a result one obtains poorer resolution, and perhaps poorer agreement with 
«3,erl»«t, in the vicinity of the shoulder of the cone than In the region 
mear the nose. The good agreement near the nose conforms with the results 
obtained In the two previous examples of flow past an obstacle. 

Figure 1-5 shows a particle plot at a late time on rtUch has been 
superimposed terschner's steady-stete bow wave position. Because of the 
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A late-tlme PAF particle plot on which has 

ti^sed ^rschn^steady-state bow f 

^ same experiment described in the caption for Fig. 





wMtaess Of the shock It Is aiffloult to pl.»olnt the EAf shock froat, 
but it seeus to agree quite well with Merschuer's, at least in the bottom 
half of the mesh. The method of Imcating the shock front in the weak 
Shock region by observing the deflection of horizontal lines of particles, 
as was done In the case of flow past a blunt cylinder, could not be used 
here. Me have not as yet been able to avoid boundary Instabilities when 
calculating flow past a cone with unstaggered particles. 

FHW THROUGH A BEMT CHANTTOT. 

The three previous examples have dealt with flow past a symmetric 
Obstacle In a straight channel, m each such case a single detached 
shock developed and moved .^stream from the obstacle, event,,.ny reaching 
a steady-state bow wave position. 

Me now consider the more complicated flow through a channel which 
contains two 90° bends. A plane shock enters the channel and interacts 
with the channel walls, producing numerous reflected shocks. The problem 

is to calculate correctly the flow through the channel while a uniform 

inflow is sustained at the mouth. 

These e;q)eriiDents were prompted by shock tube studies performed by 
He Reichenbach and his associates at the Ernst-Mach Ihstitute, Freiburg, 
Germany. In Pig. 1-6 a series of 12 photographs trace the flow through 
the Channel for one such shock tube e:5>eriment. The very con^Ucated 
structure of the flow is evident fyom these pictures which are spaced at 
40 nsec intervals. We are grateful to Dr. Reichenbach for his permission 

to use both these pictures and his detailed measurements, with which we 
compare our results. 
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Fig. 1-6, This series of photographs, spaced at 40 jisec intervals, 
traces the flow observed in a channel when a uniform 
inflow is sustained at the mouth. The pictures were ob¬ 
tained from a shock tube experiment performed at the 
Ernst-Mach Institute. 
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Fig. 1-6. (continued) 
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The experiment with which we will con5>axe was performed in the same 
manner as that illustrated in Fig. 1-6 except that a deposit plate was 
placed at the mouth of the channel (see accompanying diagram) in order to 


. . 
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obtain a controlled uniform pulse. The incoming shock had an overpressure 
of 2.29 relative to the ambient air initially contained in the 2 cm wide 
channel. This fact can be translated into the following initial data for 

the PAF piroblem: 

X component of velocity (cm/usec) 
y con5)onent of velocity (cm/psec) 
particle internal energy (10 joule) 
particle spacing in the x direction (cm) 
particle pacing in the y direction (cm) 
particle mass (gm) 
time stqp (psec) 

artificial viscosity coefficient 
(type 1-17h force) 


u = 0.0218 
1 

V, =0 

» 5.6255 
6 x^ = 0.25 
6 y^ = 0.25 
m = 1.559 X 10 
6t = 1.0 


-k 


u = 0 

o 

Vo»0 

j = 2.8155 

o 

6x = 0.555 
o 

8 y = 0.555 
•^0 


OJ = 0.1 
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where, as before, subscript 1 is shocked air and subscript o is ambient 
air. The calculations are performed in a Caitesian coordinate system. 

In Pig. 1-7 a comparison is made between the PAP pressures recorded 
at the two observing points in the channel (shown in the diagram on the 
previous page) with average curves (the dashed lines) drawn throu^ the 
experimental measurements at these two points. Dr. Reichenbach warns 
us that the relationship between his measurements and the pressure is 
not strictly linear, althou^ he points out that the deviation from lin¬ 
earity is not very great. The contparison with the PAP pressures can 
therefore be evaluated onOy from a qualitative standpoint. 

Nevertheless, the agreement between these resolutions of a very cc»n- 
plicated flow structure is remarkable. The shock front arrival times at 
observing point II (at the end of the channel) are aOmost identical and 
the magnitudes of the pressure jumps at points I and II are in good agree. 

ment. Purthermore, the major fluctuations in pressure after the initial 
jumps compare quite nicely. 

COUAPSE OF A SPHERICAL SHELL 

We do not attempt to justify the development of the PAP method on 
the basis of competition with the numerical methods already in existence. 
Several of these have been carried to a hi^ degree of sophistication 
and, when applied to problems for which they are well suited, may give 
results superior to those which could be obtained with PAP. 

Rather, the strength of the PAP technique lies in its abiHty to 
resolve problems ^ich couOd not be handled by other methods or which 
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TIME (/ts) 



Two oonrparisons between PAF pressures (mUS 11ms) ^ 

Bessurements (dashed lines) from the lOThlem 
involviiui flow throu^ a bent channel. The top plot 
S^sIS compartso^t meaMring polM 1 
and the bottom plot at measuring point II at the end of 

the channel. 






could be attempted only with great difficulty. Such problems fall into 
two categories: those involving curved or otherwise complicated bound¬ 
aries and those whose changing geometry hampers any type of mesh resolu¬ 
tion. The previous example, involving flow through a bent channel, is 
one which would be extremely difficult to solve by most other numerical 
methods because of the con5>licated boundary. If the bends had been 

curved rather than at right angles, it is doubtful that any method other 
than PAF could be applied. 

An example of the second category of problems is the calculation of 
the collapse of a thin spherical shell of material to the shape of a ball. 
This type of problem does not lend itself to easy solution for any type 
of calculational method Mhlch is tied to a mesh, whether the mesh is fixed 
in space or moves with the fluid. The primary difficulty encountered by 
the fixed cell schemes in calculating such a collpase is one of spatial 
resolution. To cover the entire region with a mesh that is fine enough 
to resolve all the details of the flow would require a. storage capability 
which exceeds that of even the largest computers. 

For the methods whose mesh moves with the fluid, this is no problem. 
Instead, the problem begins ^en the shell approaches the final stages of 
its collapse. Then, cells which were roughly square initially become so 
elongated and twisted that accurate calculations can no longer be made. 

The life of the problem can be prolonged by rezonlng distorted regions of 
the mesh as needed, but in no way can the problem be carried to completion. 
To do this would require that there be interactions between cells %diich 
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form opposite sides of the shell and, in these methods, interactions are 

restricted to cells that vere originally contiguous. 

To evaluate the ability of PAF to calculate such a collapse, a 
rather sinple test vas prepared. A spherical shell of coM particles vas 
given a velocity of uniform magnitude toward the center of the sphere. 

The material had a polytropic equation of state (with 7 = ^) so, with zero 
internal energy, only dissipative forces were in effect initially. These, 
it was hoped, would remain small throughout the early stages of the prob¬ 
lem so that the collapse could be maintained without the assistance of an 
external driving force (for which there was no provision in the code). 

The problem was calculated in cylindrical coordinates in order to 
take advantage of the cylindrical symmetry of the shell, although polar 
coordinates were used in the initial arrangement of the particles. They 
vere located along radial lines from the center of the ^here* and then 
displaced slightly from their original positions. Other details of the 

setup are listed below; 


internal radius of the shell 
particle pacing in the R direction 
particle spacing in the 0 direction 
z component of velocity 
r component of velocity 


R s 20.0 

o 

SR = 0.56 

80 0.026 « 120 

u = -cos 0 
V ss -sin 0 


%o differentiate between radial distance as measured from the center of 
iSe that measured from the axis of oyUndrioal sysssetry, « 

refer to the former as R, and to the latter as r. 
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pcu^fclcle mass 


particle internal energy 


m * product of initial 
radii, rR 

J = 0 


time step 

^ 6t = 0.1 

artificial viscosity coefficient 

(type I- 17 C force) c, = i n 


The left-hand series of pictures in Fig. 1-8 show particle configura¬ 
tions for this problem at two time unit intervals, from time zero until 
collapse. The calculation proceeded without difficulty, althou^ it was 
necessary to reduce the time step from 0.1 to 0.0125 during the last part 
of the problem when the inteiparticle spacing was small. 

The long black specks which appear in some of the last pictures 
r^resent particles whose internal energy has gone negative. They lie 
along the edges of the shell and have been cooled by the rarefactions 
^ch are working their way into the material. It is apparent that these 
particles do not attain large negative internal energies, for they do no 
damage to the calculations and in subsequent pictures are positive again. 

The pictures seem to indicate that the shell retains its spherical 
symmetry very nicely throu^out the collapse. This mai^ be confirmed by a 
study of the kinetic energy history of the problem, m spherically sym¬ 
metric shell motion the r and z conponents of kinetic energy can be ex¬ 
pressed as 


~ S J sin e dr d©\ 

0 R^ ^ ^ 
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tv - 2K y(K + 2K ) as an appro- 
One mlg^t therefore choose the function v.Ky z r z 

priate «asure of the epherlolty of the collapee; the e»re neerJy this 

expression approaches aero, the greater the sphericity. These valnes 

are Usted at too tine rmlt Intervals for this problem in the eeconJ 

colnmn of Table 1-2. Very ®x>4 sphericity is maintained nntll a time of 

12, at v*lch time the decelerating force is starting to become important. 

This leads to increased fluctuation and, hence, decreased sphericity in 


the velocity fieM. 

A logaritlmiic plot of the total internal energy of the system as a 

function of time can he seen in Pig* 1-9• 

The prohlem vas also calculated hy PAF using a force function derived 

from a "stiffened" gas equation of state [see Eqs. (l-25)]. Initial con¬ 
ditions vere exactly the same as those used in the first prohimm, p. 69, 
except for the additional specifications of normal density and the 
coefficient of the "stiffening" term, vhich vere 


p s 11*003, 
o 

a = 0*0676* 
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TABLE 1-2 


Al^ASURE OP the SPHERlEinr OP TWO lEOBIEMS 
WHICH DEPICT THE COLLAPSE OP A SPHERICAL 
SHELL OP MATERIAL 
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The normal density was chosen equal to the average initial density of 
the shell and is approximately that of silver. The coefficient, ^ch 
is the square of the sound speed under normal conditions, is one which 
vould he appropriate for a silver shell. 

Before discussing the results of this run it would, perhaps, he 
instructive to point out one limitation of this force function. It is 
possible to have attractive interparticle forces with this equation of 
state whenever the interparticle ^cing is greater than that appropriate 
for normal density material. While very small forces of this type would 
he acc^tahle for metals, they should he limited in magnitude hy the ten¬ 
sile strength of the material. No such cutoff has as yet been included 
in our calculations so that the only limit on attractive forces in this 
problem is that provided hy the maximum interparticle search radius. 

A study of the particle configurations for this problem, which are 
shown in the right-hand series of Pig. i-8, and a comparison with the 
configurations from the previous problem on the left will demonstrate 
some noticeable differences which arise as a result of this added term 
in the force function, m the plot at time 2, several particles with 
negative internal energy are observed, where none are present in the 
left-hand plot. These negative internal energies develop when the cold, 
initially staggered particles rearrange themselves in an awiempt to 
attain a more uniform interparticle spacing. Once this has been accom¬ 
plished (about time k), negative internal energies seem to be associated 
with only the rarefied regions along the edges of the shell, as they 
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vere Iti the previous prohlem# 

Another result of this early rearrangement of particles in the 
second problem is a much more rapid grorth of interparticle dissipative 
forces than vas observed in the first problem. This causes an acceler¬ 
ated growth in the internal energy of the system (as can be seen in 
Fig. 1-9), which, in turn, leads to an increase in the nondissipative 

forces. 

When one couples these increased interparticle forces with the re¬ 
sistance to compression offered by the stiffening term, it is not hard 
to understand the much greater thickening of the shell in the right- 
hand series of pictures than in the left-hand series. As a consequence 
of this thickening, the motion of the interior of the shell toward the 
center is accelerated, so that final collapse occurs about one time unit 
sooner. In the final picture in these series the shell on the right is 
beginning to expand while the one on the left has just reached final 

collapse. 

Because the shell on the right is undergoing ejpansion at time Ik, 
the sphericity of its velocity field has degenerated, as can be seen in 
Table 1-2. Up to that time its sphericity was comparative to that of 
the shell with the poiytropic equation of state, althou^ slightly 
poorer as a result of increased inteiparticle activity. 

One conclusion which can be drawn from these two collapsing shell 
examples is that such calculations can be made with PAP without appre¬ 
ciable difficulty. With the addition of an external drawing force to 
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the PAP calculations there vould seem to be no reason why truly practical 
problems of this type could not be handled. 

The second example also demonstrates that force functions other than 
those derived from a poOytropic equation of state can be applied effec¬ 
tively with the PAF method. The particular function used, however, re¬ 
quires some modification to account for the limited tensile strength of 
the material. 

F. IMPACT OF A BIUNT NOSED CYLINIiER ON A THICK PIATE 

In the examples of flow past a wedge and flow through a bent channel, 
we dealt with interactions between fluids in shocked and unshocked states. 
Satisfactory results were obtained, indicating that the PAP method encoun¬ 
ters no serious difficulty in calculating interactions between particles 
with markedly different properties. In both of these cases, however, the 
shock strength was rather moderate. We now wish to consider a test of 
the method in which the material has been subjected to a shock of infi¬ 
nite strength. 

A cylindrical projectile of cold material strikes a cold target 
plate, producing shocks which proceed into the projectile and the target 
at the same velocity relative to the interface. We wish to check the 

accuracy of the shock properties and also verify that no calculational 
difficulties arise at the interface. 

Both materials are represented by a pol^ropic equation of state 
with 7=4. The other starting data are summarized below, where sub¬ 
script 1 signifies shocked material and subscript o unshocked material. 
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z component of velocity 

u, = 1 

o 

II 

r component of velocity 

^1 ' ° 

o^ 

II 

o 

particle internal energy 

J, = 0 

o 

11 

particle spacing in the z dissection 


8z = 1 
o 

particle (pacing in the r direction 

8r^ * 1 

5r = 1 
o 

particle mass 

m = initial radius 

time step 

6t = 0.05 


artificial viscosity coefficient 
(type 1 - 17 c force) 

OJ = 1 



The initial particle configuration can be seen in the top plot in Fig. 

1 - 10 . 

At early times the two shocks produced by the impact are essentially 
one dimensional, so their properties can be calculated from simple shock 
theory. The magnitude of the shock velocity relative to the interface is 

given by 

y » ^7 (velocity of interface) 

_ 5 

- V 

The mass of material swept up by the shocks during one time unit is 
M = 2 jtp(projectile radius) v^ 

= 125. 

The specific internal energy is 

I = ^interface velocity)^. 
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PAF particle configurations from the calculation of the 
impact of a blunt cylinder on a thick plate at initial 
time (top) and at times 10, 15 , and 20. On the three 
bottom plots there have been superimposed the sh ock 
fivnt, the material interface, and the cylinder outline 
as calculated by PIC, 








so that during the one-dimensional phase internal energy is produced at 
the rate 

— = MI = 15.625. (1-29) 

dt 

We merely state, without illustration, that the early time shock 
velocities from the PAF calculation are in very good agreement with the 
above prediction, and proceed with an evaluation of late time results. 

In the two-dimensional flow regime we have no simple theory with 
which to make comparisons. Instead, we compare our results with those ob¬ 
tained using another numerical technique, the Particle-in-Cell (PIC) method 
[I-7]. Admittedly, this makes the evaluation of results scmewhat less 
definitive than would be the case if we used an actual physical e3q>erlment 
for a standard. To the ejctent that the two methods agree, however, we may 
draw confidence and in regions of disagreement we can seek explanations. 

The initial conditions for the PIC problem were the same as the PAF 
conditions, exc^t that they were expressed in terms of the PIC model. A 
cos 5 )arison between the late time geometries of the flow as produced by the 
two calculations can be seen in the bottom three plots of Fig. I-IO. 

These show PAF particle configurations at times of 10, 15, and 20, on 
which have been superlB 5 >osed the front of the shock which is proceeding 
into the target plate, the material interface, and the outline of the 
projectile as they were observed from the PIC particle plots. 

The excellent agreement with regard to the position of the shock 
front indicates that both methods are probably calculating at least the 
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gross features of the discontlimity correctly. There is some disagree¬ 
ment in the position of the interface and there is reason to believe that 
part of the blame can be laid to each method. Because of the more randwn 
nature of the PAF particle movements, one expects (and sees) a less veil 
defined interface from that method; this could probably be in5»roved with 
increased resolution. On the other hand, a detailed examination of the 
PIC quantities Indicates that there exist rather large fluctuations in 
density and internal energy in the vicinity of the interface which re¬ 
sult from the manner in which internal energy is distributed in mixed 
cells. Some of the disagreement may be caused by this difficulty. The 
outline of the projectile agrees quite well, although there is a minor 
discrepancy in the low density crater region. 

In Fig. I-n we see a comparison between the internal energy his¬ 
tories of the PIC and the PAP versions. The figure also includes a plot 
Of the one-dimensional internal energy rise rate, Eq. (1-29), with ^ch 
they both agree nicely at early times. The two methods are in excellent 
agreement up to a time of about 10, after which they diverge; the Pic 

profile attains a more gentle slope, while the PAP method continues the 
rapid rise. 

This disagreement is someehat puzzUng, in view of the fact that the 
outline of the shocked region agrees so nicely at late times. One might 
suppose that the depressed Internal energy curve of the PIC method re¬ 
sults from greater velocity fluctuations for that method, hut this con¬ 
clusion is not borne out by detailed studies of the data. A 
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of the variance of velocity magnittade for sample regions of the fluid in¬ 
dicates that the two methods eaqaerience about the same amount of velocity 
fluctuations. 

The divergence of internal energy growth rates persists even idien 
the problem is recalculated with different values for the time step or 
the viscosity coefficient or with increased resolution. This would indi¬ 
cate that the dlscr^ancy represents an inherent difference between the 
two techniques, but which is more nearly correct is difficult to say. 
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CHAPTER 5 

THE FOTURE HEVELOPMEIW OF PAF 

The development of a numerical model of a complicated physical 
system generally proceeds through several distinct stages of evolution. 
First, a computer code embracing the basic concepts of the method in the 
simplest form in which they can be expressed (usually a one-dimensional 
coordinate system, if applicable) is written and tested against uncompli¬ 
cated problems for which easily verified solutions exist. This stage of 
the development often involves a considerable amount of experimentation 
with mesh geometries, boundary conditions, and the form of the difference 

equations. 

If the results Justify it, the technique is then extended to a 
more ccmsplicated frame of reference and applied to more difficult prob¬ 
lems. Such an extension often introduces a host of new problems for the 
representation and may even require some further refinements in the basic 
formulation. Depending upon the amount of success attained in these more 
difficult problems, an evaluation can be made of the range of applicability 
for the method and the avenues in which further development mi^t be 
channeled. 
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Practical applications of the method to problems of immediate 
interest begin at this stage with the development of a working code which 
is designed for a specific class of problems. At the same time, devel¬ 
opmental work continues, extending the method to more complicated calcula- 
tioneil regime s^ 

The evolution of the PAF method has followed this pattern very 
closely. This report, in fact, marks the culmination of the second stage 
of development and we are in the process of writing a working code to be 
applied to several hydrodynamic problems for which the method seems ideally 
suited (see p. 66). Simultaneously, plans are being made to extend the 
method in order to include the calculation of the interaction between 
materials with different equations of state. 

The remainder of this chapter is concerned with a discussion of our 
futxn-e plans. 

A. THE HEtf COIffl 

The PAF developmental work has often been severely hanrpered by insuf¬ 
ficient particle resolution. For example, the problem of flow past a 
wedge, which is described in Chapter 2, was halted far short of steady 
state because there was no longer any machine storage available to accommo¬ 
date new particles entering the system. Other problems, particularly those 
involving flow past blunt obstacles, were carried to completion only by 
employing the minimum of resolution. 

The reason for this difficulty is that PAF requires a large amount 
of machine storage per particle. Not only must we store the quantities 
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needed to describe the state of each particle, such as msa, coordinate 
cosionents, the present and advanced values of the velocity components, 
and the present and advanced values of the Internal energy, hut ue must 
also keep track of the Interpertlcle relationships. For each particle, 
space 18 set aside to Ust up to four nelghhors and the Interparticle 
forces vlth these nel^hors. Also storage Is required to accommodate 
the cellular mesh mich Is used In the neighbor search routine. 

Me have attempted to minimise this prohlem hy the use of revolving 
storage vhenever possible. Several quantities vhlch are not being used 
during the same phase of the calculation and vhlch must be recomputed 
each cycle share the same memory regions. Also, the storage required 
for interpertlcle forces Is kept at the lovest possible level by storing 
only the kernel of the force function [the part vhlch Is common to F^/x), 
F (x), P (y), and FjjCy)) and by recomputing the additional factors 
whlnever ^cessary. Even vlth these economies, hovever, a minimum of 
15 memory voids per particle Is required for FAF calculations. 

Mlth this storage limitation it has become necessary to graduate 
to a calculating machine vlth a larger memory capacity in order to be 
able to apply the method accurately to a variety of problems. The nev 
PAP code is being vritten for the IBM 7050 machine vhlch has a last 
memory capacity of 98,000 vords. This added faclUty vlU permit us to 
calculate vlth as many as 5000 particles, vhlch should be sufficient 

for the problems we now have in mind. 

In addition to Increased resolution. Other improvements are 
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ejected from the nev. code. The facilities for data presentation will 
certainOy he expanded. One such addition will he streamline plots. 

These should not only he imraluahle for illustrative purposes hut also 
as a diagnostic tool. Other features which are under consideration in¬ 
clude moving pictures of particle configurations, graphs of functionals, 
such as pressure, density, and internal energy, and listings of local 
variables, averaged over small regions of space. At the time of this 

writing, several of these have already been tested and found to he ex¬ 
tremely useful. 

— REVISED FQFiM OF G?HE EMERGY EQUATION 

Another modification to he included in the new code, this time of 
a computational nature, has to do with the calculation of internal energy 
change. The equation which governs this change has been rewritten in a 
form which is designed to eliminate negative internal energies. 

In some of the examples described in Chapter 2, a few particles 
have developed negative internal energies; usually these were particles 
that were located along the edge of an unconfined fluid. These negative 
values were tolerated because the particles were either isolated, and 
therefore had little or no effect on the rest of the system, or else 
they were only slightly (and temporarily) negative, and again had 
Httle effect on the total calculation. Wo physical defense for negative 
internal energy is tenable, however, and ft-om a calculational standpoint 

it can do great harm because it leads to attractive intei^article forces 
for idiich there is no cutoff. 
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In the present code the change in internal energy which results from 
an infcerparticle encounter is shared equally hy the two particles, (in 
cylindrical coordinates this is only approximately true.) The magnitude 
of the change is therefore the same for the two particles hut the percent- 
age change may he quite different. One particle's internal energy might 
even go from a positive to a negative value while its nei^hor's was 

altered only sli^tly. 

The method would he a better model of true physical changes if the 
variation of a particle's internal energy was more strongly dependent 
upon its magnitude. The purpose of the calculational change we now de¬ 
scribe is to introduce this dependency into the energy equation. 

To accomplish this in the plane coordinate version we rewrite 


Eq. (1-6) in the form 



(1-30) 


A criterion for avoiding negative internal energies can quite easily he 
demonstrated for the case where is the interparticle force function 
associated with a polytropic equation of state. To see this we write 
Eq. ( 1 - 30 ) in difference form. 



(y - l)5t 

2 





vhere we have neglected the dissipative terms, whose only effect can he 
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is then 


to increase A sufficient condition for instiring that 8J^/j^ < i 

a j 

|ui 6t 

* * Tno xr 


■max ^_T_ 


re- 


6 r . ' f . * • 

man (7 - 1 )N 

This criterion is not violently restrictive, but is comparable to 
strictions already required for accuracy. 

Energy conservation is not hampered by this modification of the in¬ 
ternal energy equation. To demonstrate this we begin by summing Eq. 
( 1 - 30 ) over all the particles in the system. 




(i-ja) 


The sum Is, of course, ludependeut of the order of summation, so that 
may revnrite the equation 


we 


5^^ [v(“f-j)]. 


■1 


-II ■“!)]> 

0 i 


dissipative term in the force function remains 
the^nis"^? Ilf It •' difference in treatment. 
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Conibiiiing 


vhere the last expression is possible because 
this expression with Eq. (l-52), we find that 

0 ^ 

in which form the summation is now identical with that obtained from our 
original conservative equation [see Eq. (I-6)]. Conservation is assured 
in the finite difference form of the equations if, as before, the veloc¬ 
ities are evaluated at time t = (n + -g) St, 

A similar treatment of the cylindrical coordinate equations does not, 

however, produce energy conservation. In that coordinate system the analo¬ 
gous equation for the change in internal energy would be 


dJ 

—1 = 
dt i 


,« 1 J 


1-55) 


When we sum this change over all of the particles in the system and add 
to it the total change in kinetic energy, we find, after considerable 
simplification, that the total energy change fails to equal zero by an 


amount 

i = 51 Z* [C-^) • 

i i 

To correct for this discrepancy and insure conservation, we therefore re¬ 
write Eq. (1-35) in the foim 
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Once egein, this conserretlon carries through to the difference equations 
if the velocities are time-centered. 

At the tine that this report uas urltten several tests of the neu 
code had been made and these prelliDlnary results Indicated that the re¬ 
vised Internal energy calculation tos perfomlng satisfactorily. In com¬ 
parison ulth runs perfoimed vlth the earUer version of the code. It iras 
apparent that the large scale features of the calculation vere essentially 
unaltered, but the absence of negative Internal energies uhere they had 
been present before testified to the effectiveness of the neu technique. 

It remains to be seen whether this revised calculation will obvUte the 
particle destruction procedure which liad been necessary In order to avoid 

boundary Instabilities In certain cyllndrlcally syetrle flow problems 
(see pp. 42). 


c. MUIIFI-MATERIAL faf 

The new PAP code is a one-material program, as is the one described 
in this report. Plans exist, however, to extend the method to include 
interactions between materials with different equations of state. 

Fundamental to these plans Is the concept of a force function between 
particles of unlllce kinds. In the absence of a correspondence theory 
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(such as that descrihed in Eef. 1-2) for this type of interaction, ve de¬ 
rive the interpaxticle force function on the basis of preserving fluid 
equilibrium. In particular ue require that the average pressure be con¬ 
stant on either side of a material Interface in a fluid in equilibrium. 

1. Plane Coordinates # • • X X X 

For simplicity, let us as¬ 
sume that the psirticles which ® ® 

represent this fluid are arranged ^ ^ # I X X X 

in a square array (see acccmtpany- | 

INTERFACE 

ing disigrsiJ®) sincL tlis/b the cooirui"* 

nate system is Cartesian. The derived conditions should then be necessary, 
but perhaps not sufficient, for general equilibrium in a plane coordinate 

system. 

The material to tte left of the Interface we shall call dot material 
and that on the right x material. The Interparticle spacing is the same 
for both materials; hence, the dot and x particles must have different 
masses if their densities are to be different. likewise there are two 
distinct interparticle force functions, 




for the two 


’ij x' iy i' y 

fluids. At the interface, equilibrium requires that these 


two be equal: 



If the two fluids 


were polytroplp, this equation would be written 


- '>J. - OJj, 

For we try an expression of the form 



and solve for a value of 1 which is consistent with the equation above 





J + (• J J 

1 ) J = X 


• 

2 


2(7^ - i)(r^ - 1 ) 

7 T~y ITa 

• X 


The interparticle force function between two fluids with polytioplc equa 
tions of state can therefore be written 




r o-t - I)(r, - 1)(J, w,) 

L in t rj - 2)r^j 


notice that when 1 and J are the same material this expression reduces to 
the usual one-mterial force function. 

If the two fluids had stiffened gas equations of state we could write 


Eqs e (I- 35 ) as 


(y - l)J m a 

' ' ^ mm 


(y _ i)j m a 
X ^ X X 


Using the same value of X as vas used for the polytropic equation of 
state, ve can reduce this equation to 

The only symmetric value for T is the average of the two terms 


m a + m a 

. . XX 


“ m a + m a J 

*- , , X X 


The force 


function appropriate for two unlike stiffened gases can there- 


fore he written as 


(y. - iKr, - i)(Ji + J,)' 


' =1J 


/ 

,1 L (ri“ 






Vi * ”j'‘j 


2 , Cylindrical Coordinates 

Consider the problem of maintaining equilibrium among a square array 
of particles in a cylindrical coordinate calculational system. If we 
confine our attention to the z component of the inte 3 rpartlcle 
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relationships, we see that the problen is entirely unchanged froai the 
plane coordinate case. Therefore, the plane coordinate form of the 
interparticle force function will apply as well to the s component In 
cylindrical coordinates. 


The r component of this interparticle force requires additional con¬ 
sideration, however, because of its asymmetry. We begin by recalling 
that the one-material form of this force was derived subject to the con¬ 
straint that an unconfined fluid of uniform density and internal energy 

must have the same rate of ej^ansion in the positive and negative radial 
directions. As a consequence of this demand, a confined fluid in equili¬ 
brium which is subject to the proper boundary conditions must remain 
forever in equilibrium. We wish now to extend this requirement to the 


case where more than one fluid is so confined. For simplicity of calcula- 
tion we again specify that the particles be arranged in a square array. 

Consider a radial line of these ^ ^ 

particles (see accompanying diagram), 

X, 

numbered 1, 2, 5, and 4, which crosses? _ ^ _ 

INTERFACE 

this interface. Particles 1 and 2 are ® g 

dot material, 3 and 4 are x material. 

In the one-material case, when the den- ^ 


sity and specific internal energy were uniform, it turned out that the 
force on a particle from its neighbor above was exactly balanced by the 
force exerted by its neighbor below. We now require that this hold in 
the vicinity of an interface, so that 
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(1-36) 




The one 


-material analogy is carried even further with the requirement that 


TP __2 F . 

^25 rg 32 


(1-37) 


If a force function can he discovered which satisfies conditions (I-36) 
and (1-37), then equilihrium will he assured, for the pressure will he 

uniform in the vicinity of the interface. 

If the two fluids have polytropic equations of state, condition 

(1-36) can he expressed 

^25 6rir + rj^) Br 


^ (7 - 

-El • C. • • 

F32 -- 5 r 


where we have made use of the fact that in a uniform density fluid the 
mass is proportional to the radius. We now use condition (1-37) to deter- 
mine the form of the interparticle force function; 


F 




25 


6 r 


6 r 


J 2 

'3Lbr(r, + r Jj 


(.7. - iXr,- 1)(r2 + Tj) 


- 9 T- 


< u 



Generalizing, now, to any array of particles we write the interparticle 
force function as 


r - ’>(7, - i)r -i(j,+ 

■ ^17“ - o J^ 

(Z) . " Ji) 

L J ry • 


The uniformity of pressure along the radial line may now be veri 

fled by applying either of the relations from conditions (I- 36 ). For 
example. 



' ( 7 , - i)(y^ - i K ] (J, + J,) 


r 



1 )r_m I 
3 X X 


8 r 


r = 



simplifies to 


iy - l)m I = ( 7 ^ - l)m I , 

so the same pressure exists on either side of the interface. 

Applying the same technique to two fluids with "stiffened" gas equa¬ 
tions of state in cylindrical coordinates, one can obtain an Interparticle 
f*orc© function with coniponcnts 




- J 

1 

r- 1 

1 

•rl 



.L J' 



■’■ L J L ■*■ 


•" Voi\^ 




]}• 


It can be verified that, vhen particles i and j have the same equation of 
state, these expressions are identical with Eqs. (I- 25 ). 

Notice that the kernels of these cylindrical coordinate multi- 
material force functions are not identical for the r and z components, as 
they were in the single material functions. As a result, the multi¬ 
material PAF code will require twice as much memory for the storage of 
interparticle forces as does the present code. 
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apeendix a 


V2RIAL APPROACH TO CORRESPOMDENCE THEORY 


Neglecting the dissipative forces, we may write for the PAF particle 
dynamics the eqiiations 


du 


m. 


5t " I 


dx. 


d-r = 


Prom these can he derived 

^ (m.x ‘u'.) - m u ^ •t f 

i 

Suppose that the particles are confined within a rigid circular (or spher¬ 
ical) region hy means of images beyond a wall of radius R, The above 
equation can be summed over the particles within R and rearranged 


as 


follows: 


I ® -1=I X 


J,int 




J,int 


,},all 1 


■I Z' 

J,ext 1 


x.*B. .f... 
J ij ij 
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The sum over means the sum over particles within R, while the 

sum over J,ext means the sum over the exterior images only. Summing over 
J^all includes both. 

Consider now liie long time average of this equation. The first term 
vanishes, since m.S^.'u. is hovinded for these confined particles. The 
second term vanishes because the PAF correspondence is supposed to force 
all particle kinetic energy to be macroscopic only. Thus, we obtain 

I Z E ' Z Z ’ 

J,all i J^ext i 


where < > means the appropriate long time average. For the exterior 
particles, ^ vh&re is the distance out to the image from 

the wall. In n dimensions, is the surface pressure, where s is 

the mean interparticle spacing. Thus, with the assun 5 >tion of strong corre¬ 
lation, the virial equation becomes 

1 (N^N* + N^)sf(s) = Ng(R + ' 

vhere is the number of interior particles, is the number of exterior 
images, N* is the number of neighbors per particle, and p(mp) is the pres¬ 
sure as a function of mass density. It can be shown that ^j^s, so 

that 
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Thus, in order to have a correspondence equation which is independent of 
necessary that N = 2n, proving, therefore, the contention 
that the numher of neighbors per particle must be twice the number of 
^ace dimensions. We also obtain 

which is the same as the result obtained f3x)m the liouviUe correspon¬ 
dence theory — see, for example, Eq. 
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APPENDIX B 


THEORY OF CORRESPONDENCE FOR CYLINDRICAL COOEH)INAIES 

We vish to extend the correspondence theory to show that three- 
dimensional problems with cylindrical symmetry can be approximated by 
a two-dimensional PAF model. In other words, we ^>rLsh to use a two- 
dimensional array of particles and forces between them to represent a 
three-dimensional array of particles and corresponding forces. 

One of the restrictions we should like to make is that the force 
between two particles be dependent only on the plane coordinates of the 
particles and the usual variables 9 and i?. (Notation throu^out this 

appendix is the same as that in Ref, 1-2.) 

One possible scheme suggests itself immediately. V7e could con¬ 
sider each two-dimensional particle as the representative of a ring of 
particles. Then the force between two such representatives would, be just 
the average of the actual force over the entire ring. Unfortunately such 
an association of the two-dimensional array of particles with rings also 
means that a "self-force” would be imposed on each member of the array 
because of the other particles in the ring it represents. This 
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self-force leads to considerable trouble, as »111 be seen below. However, 
if the Initial distribution were such that this self.force could be 
ignored, this approach is partially successful. 

A slight alteration of this technique is to assume that the initial 
distribution is such that, while it is not formed of rings, little error 
IS involved in computing interparticle forces by means of an average. 

However, it seems difficult to justify this procedure on anything other 
than an intuitive basis. 

It should be noted that it is necessary to assume that variables 
such as p,, ^ and u do not change with e (where r, e, and t are cylindrical 
coordinates and the t ails is the axis of sjmmetry). These assusi,tlons 

are justified by cylindrical symmetory. 

We begin with N particles in a three-dimensional array whose motion 
is governed by a distribution function where j runs from 

1 to N and r^ is the ordinary three-dimensional position vector. 

We first integrate out for each j to obtain a new distribution 
function 

0 

where is the ordinary two-dimensional position vector. 

This means, of course, that we are considering each particle only 
in terms of its r and z coordinates. As in the previous correspondence 
[1-2], we conclude that 
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I ^ ■■■ 

V 

for a generalized volume V whose houndary is not crossed hy particle tra. 
jectories. This leads to the IdouviUe equation 





(B-2) 


(I - that the particles are intrinsically identical (l.e., have no 

Identifying property such as mss vhlch they carry along In space), then 
the reduction shown on p. 15 of Bef. 1-2 is valid and the result Is 






I) + heat flux terms] 
(B-5) 

X r dRg dUg dcpg. 


As before, ve assume that 


c(?pe(?2 - J'l-Pip.'Pu) ®(“2 - - '^12> ^^'^2 - “"l ■ 
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Since = 0, this implies that = 2nz^ and 

1 r 2jt o 2jt 

0 0 


= pCr^) 6 ( 5*2 - ^ - “12) 6(CP2 - 9, - 
r 2:r 

^ ^ ^ ^0 ''^"’l2^®2rPl2^'Pl2) *^^2 ■" ^^2 • 

Note that ‘^(ri2"Pl2"'^12^ depends only on 0^ - 6^ since 


(B-4) 



^•,2= (^2 - 

Since 

and are : 

r2jt 

r 2jr 

L - 

/ Odd, de3 = 


Since 


+ r^ + - 2r^r^ cos(e^ - e^). 


■2jt r p2n-d^ 


<*( 9, - fg) a{ e, - e^)] ae^ 


r 23r 

= 2rt / o(e) d(e) s 27ta' 
O 0 


so that 


o »2 = 4) 6(?2 - ?, - ?, 2 ) 6(93 - 9 , - 6 , 3 )., p.(?3 - VP, 2 , 9 ,,) . .. 


Substituting Eq. {B-5) Into Eg. (b-j) and Integrating vlth 
respect to du.j dcp^ gives 

r i ■* St *- ^-.<P5) = 0 , 


(B-5) 


(B-6) 
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where 



This is just the equation of continuity in cylindrical cooralnates. ^ 

In fashion, integration ulth respect to ? d? dcp and <p dSdtp 

yields the equations 

du A V (B“T) 

p ^ + p(u • v_pu = pP, 

dt R 


A (p5) + V_>.(pu u) +p^= pP 

^ R 

and 

^ + V_, ) + P I 

R 

which are the ej^tected momentum and entropy equations. 

j. ,.^13 ACilA-¥ — iT = 0 in order to derive 

Note that one need not assume R d6/dt - ^ 

these equations, as long as ue have (S/^d) u^ - 0. 

The ahove is a stral^tforuard adaptation of the previous theory. 

Me must no. choose a fom for the force function between two particles, 

and then evalmte the expressions pP and pQ. 

m actual practice, the twcdimenslonal array of particles have 
mass proportioned to their initial distance from the axis of symmetry. 
For purposes of this corresEOndence, it seems much easier to assume that 
all have the same mass (otherwise, the basic assumption that particles 
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are identical would be invalid). 


As noted above, if we assume that each 


particle is the representative of a ring of particles, a self-force «st 
be taken Into consideration. If one goes back to Eq. (B-j) and lets P^, 
contain ? elf i"aapendent of particle 2 , then It can be readily^' 
Checked that the momentum equation will become 



= pP + pf 


self* 


(B-9) 


y symmetry this fgg^f ^ radial force and will in general Tna i c e it im^ 

possible to express the right-hand side of Eq. (b- 9) as the gradient of 

a function. 

In addition, this Is nost reasonably expressed as an average 

over the hjq^thetlcal ring. In other «rds. If is the usual force 
between particles, one expects the self-force to be something like 


_L r2jt ^ 

Jo 21 


and this may well diverge. 

Therefore, we must discard the self-force concept. However, we 
wish to preserve the averaging technique. To do this, we appeal to 
cylindrical symmetry to argue that the averaging procedure Is the best 
two-dimensional description of the actual three-dimensional force dis¬ 
tribution. Hence, while each particle is not now regarded as the rep¬ 
resentative of a ring of particles. Its e coordinate Is essentially 
arbitrary. 
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Therefore, we 


talce as the force function "between two particles in 


our two-dimensional array the quantity 


1 P 

^21 ^^2 
0 

’■-[h Jj" ^ <*2 ®2 • ‘®} 

_ fjL r f(r ,cp,») (z - z ) delz 

l2ii Jq ^ 12^^12^ ^12 ^ ^ -* 


' ^21,r^*^2t,z"’ 


where f^^ is the usual inteiparticle function "between three-dimensional 
particles (for convenience, only nondissipative forces are discussed 
here — dissipative forces will "be discussed below). Therefore, 


= B ^ 21 , z >^2 '^2 


(B-10) 


= B //p(B2’‘'’^21,/2 ^2 


(B-n) 


One now makes the usual first order approximations and applies symmetry 
considerations to ottain 

( f^>z = - ^ //‘’■ f21 , z<"2 - ^’“"2 '^2 ® J ' 


[^//‘’’^ 21 , r '’^2 - =^ 1>‘"2 ^=^2 ^^ 2 . • 


(B-15) 
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This would be fine if the two quantities in brackets were in fact 
equal. The bracketed quantity in Eq. (B-12) can be e3q)ressed 


J 


12 ^ 


where dr = r^ dr^ dz^ de^ is the three-dimensional volume element. How¬ 


ever its counterpart in Eq. (B-13) is 


^ 1 °' ^ - ®,) (^2 - r,) dT. 

In general, these expressions «111 not be equal, since o' can be a com- 
plicated expression in r^ and z^. 

The problem is that we have lost a good deal of symmetry in these 
e;q)ressions by integrating a and f separately over 0. Initially, a is a 
cutoff function symmetric in all three directions, which allows one to 
make general first-order e:q)ansions. However, a* is a nonsymmetric cut¬ 
off function in two directions which allows one to make first-order ex¬ 
pansions only in r and z. 

Hence, this simple averaging technique for interparticle forces Is 

not valid either. Instead we must take a weighted average In order to 

Introduce sjnnmetry Into the appropriate Integrals. Therefore, we now 
define 


f = — '* 

21 2it 

f o de^ 


n2n 


-no- 


so 


that Eqs. (B-10) and are replaced By 


and 


' -»/ ‘’^"2’ ^ ‘"2 ■*"• 


(B-14) 


(B- 15 ) 


Now we can make three-dimensional expansions, using the fact that 
0 is a symmetric cutoff function, and the result will Be 

Which are equal By symmetry. Hence we make the association 

P=4/ 

.here p is the pressure. Exactly the sa« syneietiy consideratlcns enter 
vlth regard to dissipative forces. If usual dissipative f 

between three-dieensional particles, one takes as the dissipative force 

between particles in the two-dimensional array 

r2rt 

j Sg2, IV 


i _ 


21 2n 


J o de. 
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and the resulting Integrals pP and pQ will he analogous to those evaluated 
in the previous corre^ndence theory and in Appendixes C and D. Hence, 
one will get derivatives of a quantity interpreted as the stress tensor in 
pP and the heat conduction term V»(/c^) in pQ, 

To conclude, it has been shown that correspondence can he achieved 
for problems with cyHndrical symmetry when all masses are assumed equal. 
Attempts were made to relate this ^proach to that actually used in calcu¬ 
lations by coalescing particles in order to approximate a distribution 
with mass proportional to distance from the axis of symmetry. Thus far, 
such attempts have not been successftil. 
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APPENDIX C 


VISCOUS COERESPONDENCE FOR PAP 


As on p. 20 of Ref. 1-2, ve vlsh to find a formula for dissipative 
forces between PAF particles for which angular integration gives us a 
term to be associated with the derivative tensor 



If p(?) is constant in space, then Eq. {2k) of Ref. 1-2 is in such a 
form. However, in general, there is no reason for p to be constant; in 
fact, in almost all problems studied, it can vary considerably. There¬ 
fore, it is necessary to look for new terms to put in Eq. (22), Ref. 1-2 
Actually, several terms were omitted in the derivation of Eq. (24) 
of Ref. 1-2, and, before looking for new forces, we will discover the 
correct e 3 q)ression for the forces used in Eq. (22) 

The expression 


can be e:q)anded to 
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/[p(S;) ^ (?2 - {<’(p(?,),r,2l * I (?2 - 





?,0) 


+ ^“,a-v )" ^,(r,34,,5,o)l a?,, 

12 

where use has been made of the fact that ?and b^/d<p are e^cpected to 
vanish for a = o. We have also assumed that o was independent of cp. 

Removal of this assumption changes none of the results presented here, 
but merely adds complication. 

If a is e:5)anded to second order, we get one first-order term and 
five second-order terms. In the event that ^is linear in a [so that 
the term in (oc V) g vanishes] and that the first-order term vanishes, 

there will be four second-order terms to account for rather than one, as 
is shown on p. 20 of Ref. 1-2. 

With ?. t the first-order 

term does indeed vanisl^ and (a«V^^= 

First consider 

^21 " - ?,).V]?+ 1 t(?2 - 

the following are all second-order terms: 
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( 1 ) 




= 1 / 


^2- 

ffAr^x? — 4^ sin e dr de dcp 


MU 


4 "i 
iPoT ar I ^ 


(Here, as in the foUowing, i denotes the vector Index (no smn.ation) 
and sums are always taken over indexes k, i, m, etc., whenever re 
peated. Xj "■eans - ?,)j, and all derivatives of are at ?,.] 

( 2 ) 


' I is ^ 
-f /*is4^sr 


_ ^ 
" 5 


ha ^ ^ 1 

*SS“’ 


= ^(/” 
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(4) 


2 J 




= iSB. 




Hence, using these terms in the ejipression for p? gives 


with 


= s£ r 

i 


4 

A0r dr. 


which is not of the desired form* 

Next consider the term 

B(ri2,q),2)?2i(^r?21>- 

Again there are foxir second-order terms, and the final result 

2 a f^± ^ 


Cj “ik 


ha (^± H ^ N 

+ p * »ik sp* 

^ r 2 /^i ^ ^A 

' 3^ L” '* w * ST * ®it sr) 

i I 


where 


■hir 


4 ^ 

oBr dr 
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Since this means on]y one of the independent coefficients has 
been determined, the problem is to find another es^ression to be added 
■^Q ^ vhose integral can also be associated with the derivative of a 

•■4 

stress tensor. We hoped to find another expression linear in a; however, 
none was found to be suitable. 

Hence, attention was turned to other e^gnessions in the values 
u(i^). We first observe that 

— — ^ \ 

. +5 t (?2 - 

to first order. If we suppose that 

were part of the force, then it can be checked that the usual inte¬ 
gration and multiplication would yield an additional term 

with 

'i i(/J • 

Similarly, the term 
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(H (?,)? .[!(?,-?,).v]|^^)}. 




yields the ejjpression 




[••• 4 * 5 )]. 


where 




We note that both of these terms are of the desired form, and if g 
were equal to their sum, we would have 


(o^)j ' ^ [“S ^ ®ik 


as desired. 


Hoover, neither (V?) (^4) nor ^ can {in general) 


i 

=»/ 


be approximated to first order by the values u (r^). 

Hence, we abandon the restriction that dissipative forces be of a 
two-particle nature and look for analogous e3q)ressions involving a par¬ 
ticle and its neighbors, in the hopes that these more general e 3 q)res 3 ions 
can be approximated by ?(r^). 

The natural candidates are, of course, esgpressions in 


(vu) (-!_*UL) 
^ N + 1 ' 
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(where * are the positions of the N nei^hors of r^) and 

2 N +1 


^ ^ N* + ^ ^ 


In both cases, we are interested to see if derivatives of the vector 
field u at the midpoint of H* +1 particles can he approximated hy the 

values u(r^) to first order. 

We first observe that approximations of the form 


5(?i) + 3t?,) ' 2 ? (^ 4 % 


are only good to first order. Hence, the natural procedure of using 
such estimates for u at hypothetical middle points is not valid, since 
the result can only he a zeroth-order approximation for a derivative. 
However, the approximation 


is good to first order in the derivative, so we can assume that such 
directional derivatives at the midpoint of any line segment are known. 
Also to first order in derivatives of ii we can approximate 



therefore, if we can find equations giving us the quantities 



from the known directional derivatives at the points —i, then we will 

indeed be able to approximate derivatives at the midpoint. In one dlaien- 
sion this is trivial, since all derivatives are in the same direction. 

In two dimensions, first consider four particles at r . ? ? and 

r^. Unknown are the 12 quantities: 




¥ \—r^J ^ ^i j 

for i^j s= i < j. If we dexiote 

5? (?j) - ? (?j) by (= 

and 

by b^^. 

then we have the six equations 

Also both and satisfy the equations 


^12 ■*■ ^34 “ ^14 ^25' 


‘‘la ^34 “ ^3 ^4' 
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But this seems to exhaust our toouledge of the system end ve have only 
10 equations In 12 unlmovnB, so of course there is in general no 

solution* 

W€ now go on to five particles in two dimensions, using the 
same notation as above. It can he checked that both and must 

satisfy the five independent equations 


^12 ^54 “ ^15 ^24 

= Zl4 + =^25, 

^ 15 ''" *23 ' ^15 ^15 

■ '12 * ' 35 ' 

and 

^3 ^45 ' ^15 ■" ^54. 

These 10 equations together with the 10 eqjiations 

®id “ ^i/i,3 

allow one to solve for the 20 unknowns x^j^^ and and thus to obtai 
desired derivatives. This syafcem of equations can he easily simplified 
to the prohlem of inverting a 5 X 5 matrix, hut heyond this point further 

simplification seems difficult. 

The analysis is quite similar in three dimensions. Here we write 
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K = a. .X. 
ij ij 


+ b y 

i/ij 


c z . 

IJ ij 


As with -four particles in two dimensions, the use of six (or fewer) 
particles in three dimensions does not provide enough information for a 
general solution. However, with seven particles, the x , and ? 

ij^ ^ij» 

must all satisfy l 4 independent equations of the type 
^12 = ^13 \k 

This gives 42 equations which, together with the 21 equations K . = a.. 

^i/ij ‘^ij^ij’ Sive 63 equations in the 63 unknowns. Solution of 
this system of equations can easily be reduced to the problem of inverting 
a 21 X 21 matrix, but further reduction seems difficult. Hence, in prin¬ 
ciple, one may always approximate (to first order) derivations of a vector 
field at the midpoint of N + 1 particles only from knowledge of the field 

itself at these points as long as N* > 2n (n is the dimension under con¬ 
sideration). 

In practice, one would never invert such complicated matrixes as 
those above, but just use n sets of two particles whose midpoint was 
nearest the general midpoint, approximate directional derivatives at the 
general midpoint by directional derivatives at such two-particle midpoints, 
and use the resulting n x n matrixes to calculate the fixed derivatives. 
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Va Buat nw aho* that the gaaeralited t force above can indeea be 
fitted into a correspondence between fluid dynarfca and PAP. We retreat 
to the generaUted Ilouvllle equation on p. 12 of Ref. 1-2 



=-(y 

dt B ViLj 


^3^13 * ®3 




Where t 1«!«“4s also on the N* neighbors of the 3 particle, hut has no 
sinple resolution into two-particle form. Integration as on p. 13 of 

Eef. 1-2 gives 


^ t Sj-v^s, t j'd?^ ( 


® 21 ^ 21 ^. 




[/ 


clx_***clr* d\?*‘*du* -t-g 

2 /+i N +1 N +1 


dcp2***d9 * 
dcpa 


z * - 
N +1 

m 




] 


1 ^ fz ^Clr**‘d(p =0 

+ (n - 'ill ^ J N dt 2 N 


-me first three expressions are unchanged and the last is irrelevant for 
t.e present discussion, so it is only the fourth term vith which we con¬ 
cern ourselves, i.e., we must suggest a form for z^*^^ vhich upon inser¬ 
tion of mean values for and cp allows us to perform angular interactions 

on the term in bracketse 

we have only to prescribe a fonE for in the region where % 

Is nonzero and, in fact, we assume 
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' ^ n't! ’ N*+1 ’ ’^12'"'''l(ir*+l)J 

X p(?3,«-p(? t)1f’ . (?^. ?,.V)S;)} 


X B[cPt - (^1 + ^ - ?1*V9^)] 


in this region. As before 


= p 6 (x? - 6(9 - 9 ). 


If P» is the term in brackets above, this means that 


* 

1 P N +1 rp 


%*+1 


+I 


(N*+1)J ®1 


With the usual first-order expansions for p, a, and ^, where 

k#s 2 

* b»2 

k^2 N +1 ^ 
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it is straightforward that 


where 


and 


This is 



an e:^ession of the form ve wesre looking for. 
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APEENBIX D 


HEAT TRANSFER TERMS IN THE RAF MODEL 


It is desired that heat conduction be Incoiporated Into the PAF 
model. This is done by allowing for heat flux as well as work flux terms 
in Eq^» (4) of Ref. 1-2. Here we now write 


1 V* 1 r’* v'* 

dt " ■ 2 I ^ij -dt^ 2 L 

where H will be specified later. Among other requirements H must, of 
covirse, satisfy h(J^,Jj) = -h(Jj,J^), so it is indeed a flux between 


particles. 


In the previous discussion, the variable (? has usually been left 
unspecified. We now identify (p with internal energy. Therefore, the 
last expression in Eq, ( 6 ), Ref, 1-2, becomes 


^ ^21 ^ * 4 -(“2 - ^1 ) * « Ja'-'i )] } 


and so instead of Eq. ( 15 ), Ref. I-2,we write 
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For the term we calculate separately for all the e3i5>ressions con¬ 
sidered as candidates for : 
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(2) B(r^2,9 + I 


s‘(/;.:-)[^(^.%.vSi)]- 


(3) 0(r,2,5 + 1 B,^) V.? (5-^) gives 


4 

Car dr 


•)sf5- ^(/ & (v s^) 


DC^ig^tp + ~ 


PjP {(5^ - + [(r^ - ?,)-7] 


(/o ^ • 


gives 


(5) Similar results would hold if the generalized "g force, with terms in 


^ N + 1 ^ 


were used, since the calculations sure ajisilogous. 

Hence, if we use a combination of these forces for such that 

then what we have here is Just the e:^ession 
m ^ ^ 

. — P y.. + —J^ 0 . 

p dXj Jk 

Prom p. 123, Eq. (l4) of Ref. 1-11, 
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dt Sxj \ Sxj/ ’ 


vbere I la apeclflc Internal energy. From ®i. (19) of Kef. 1-2, 
^ + V.(^j) = P ^ * "fij ^ *^2 

since 


PiJ “ °ii3‘ 

Therefore, as J * ml, we can associate the easpression 

£ r dr- poH 

m J 2 

with 

^(x^).v.(gv?). 

m fact, since temperature is a function of internal energy, we take 


SO 


+ irl(?2 T(?,). 

that after the usual calculation we find that 


I / d?2P0H(J2,J,) = 


sr 


[4 


and we can take 


(/o 


^ r r 4 

= 5m 
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It has also been suggested that the variable (p be iirteii)reted as en¬ 
tropy, so it is of interest to note that heat conduction terms can also be 
Incorporated into PAF under this interpretation. Here we add a flux term 
to Eq. ( 5 ), Eef. 1-2, so that 



Equation (49.4) in Ref. i-i4 gives 

(ll ^ ^ “ik . 

Since ve have 



these equations are in agreement. Unfortunately use of this procedure 
requires knowledge of the temperature as a function of position. Excqrt 
for a polytropic gas in which T « I, this may be unknown. 
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ahendix e 

MASS EEABRANGEMEUT IN CYLINIDiRIGAL COORDIHAaCES 


The initial particle configuration in cylindrical coordinates is 
generally arranged in such a vay that particle mass is proportional to 
radius throu^ut a region of constant mass density. Subsequent circula¬ 
tion of particles through the mesh, however, can convect massive particles 
towards the axis and li^it ones away from it, with a resulting loss of 
resolution or creation of instability, respectively. One effective method 
for relieving the instability difficulty, described in the text, involves 
destruction of low-mass particles which have moved far from the axis. 
Another method would do just the opposite; through a mass rearrangement, 
the mass of any particle could be changed every cycle, so as to continue 
the proportionality of mass to radius. Thus, for example, we could add 
to the procedure the interparticle mass flux calculation, whereby 



where, for mass conservation, we would require "^ji* ^ leasonabl 
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choice for K . seems to he 

i j 



and it can he shown that a necessary condition for stability is that 

0 < I < 1/2. 


Now it is necessary that such an interpaorticle mass flux ix>t result 
in any net macroscopic flux of mass, momentum, or energy. This means, 
then, that the particle coordinates must change in order that subset cen- 
ters of mass remain fixed. Ifc can he shown that the only way to accom¬ 
plish this is through the simultaneous solution of the following for the 
new coordinates, r^*, in terms of those before mass redistribution, ? : 

Zj Y-i*P * 

i 



where the p modifies the sum to include only those terms for which > 0. 
(This equation is derived by insisting that all the mass associated with 
a given particle have the same center of mass adCter the redistribution.) 
This inhomogeneous summation equation can be solved by successive substl- 
tution of the first term on the right into the second. In many cases, a 
single substitution probably will suffice. 


Finally, we must redistribute energy and momentum, and the appro 


prlate equations are 


•7'^i ‘ "Aj ** * ("j • 

”r®j ’ Z*^ * (’’3 ■ 
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PART II 

SOME BASIC PROPERTIES OP PARTICLE DYNAMICS 

introduction 

The results of a pioneering study by Pasta and Ulam [II-i] indicated 
that it shou:id be possible to develop a useful technique for high-speed 
computer solutions of problems of fluid dynamics in yAUch fluid elements 
are replaced by discrete interacting particles. Their work did not, how¬ 
ever, clarify the nature of the correspondence between the model and a 
true fluid. There was, for example, no indication of how the interac¬ 
tions should be chosen in order that the model possess the equation-of- 
state properties of the fluid, and in addition there was no provision 
made for including the dissipation necessary to remove reversibility, 
later work [II- 2 , II- 3 ], however, has demonstrated that these objec¬ 
tions can be overcome. Thus it has been shown by a statistical analysis 
that it is possible to establish the relationship between the inteiparticle 
force function and an equation of state for the fluid. Also, it was 
shown that the addition of an internal variable for each particle and 
inclusion of a dissipative term in the force function can satisfactorily 
remove the reversibility. Part I of this report shows the most recent 
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status of the method resulting from these generalizations. 

We have found, nevertheless, that considerable information vhich is 
pertinent to such a generalized model can be obtained from some simple 
additional numerical experiments using the original Pasta-Ulam procedure. 
The results of these experiments are of considerable aid in the interpre¬ 
tation of the earlier work and also show how to resolve some ambiguities 
in the development of a workable generalized model for solving nonsteady 
fluid dynamics problems in several space dimensions, using hig^ speed 
computers. It is, therefore, the primary purpose of this Part to present 

these resvflts and their Interpretations. 

An additional result of these simple numerical experiments has been 
to demonstrate a somewhat different type of problem to which this sort of 
calculation can be successfully applied. Apparently it is possible with 
very few particles to be able to make meaningful statistical analyses of 
the properties of true nwlecular assemblages with a high degree of 
accuracy. Some preliminary discussions of this interpretation are also 
included, and we shall soon present elsewhere the detailed results of an 
extensive investigation. 


- 137 - 



CHAPTER 1 


THE NUMERICAL EXPERIMENTS 


We have examined in detail the dynamics of a one-dimensional set 
of N particles, interacting in pairs, and confined to move between rigid 
walls (accomplished by holding the first and last particles identically 
at rest). As in the Pasta-Ulam calculations, the force between particles 
was repulsive and inverseOy proportional to their spatial separation. 

The initial conditions for all examples were such that the particles were 
uniformly spaced with separation s, and all but the end ones were moving 
to tne right (positive x direction) with the same velocity, U. 

The equations which govern all details of subsequent motion are 



dx, 

5t^ (II-2) 

These describe the time rates of changes of velocity, u^, and position, 

x^, of particle The constant parameter, c^, is the ratio of the 

force constant to the particle mass, the latter being the same for all 
particles. 
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To derive ad energy Integral of the differential equations, multiply 
Eq. (Il-t) by u and SUB over all particles. After some manipulation, the 
result becomes 

N N-1 .v - X \ 

B-M 

Use lias been made of the fact that particles #1 and # are held identi¬ 
cally at rest, so that and axe constants. The first sum represents 
kinetic energy, while the second is potential energy. In the latter, the 
incorporation of s shows a choice of integration constant malcing the 
potential energy vanish when the particles are in their initial spacing. 

E is the (constant) total energy. (All energies are specific, that is, 
measured per unit mass.) The energy considerations are of importance 
for several purposes, particularly for judging accuracy of the computer 
solutions and for deriving the fluid-dynamic properties of the particle 

system. 

For the computer solution of the particle dynamics, the time deriv¬ 
atives are replaced hy finite difference approximations: 


du. 

_ ^ 

dt 


u ^-^1 - u ^ 

J_ a 

6 t 


dt 


n+1 n 
6 t 


in which the 


index n counts cycles of elapsed time, each of duration 5t. 
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Thus, Eqs. (ll-l) and (lI-2) become 


n+1 

u. 

J 


n 2 
u. + c 
J 






(II-4) 


n+1 


n 

= X. + u 


n+1 


5t* 


(II-5) 


In this way, then, evolution of the configuration throu^ a time cycle 
is computed by substitution of available data into Egs. (lI-4) and (II-5) 
for each particle. Note that the use of the newOy computed velocity in 
Eq. (II-5) implies that substitution into Eq. (li-4) comes first; the 
reason for this particular ordering is that it contributes significantly 
to the difference equation stability. (See, for example. Ref. II-2, 

PP. 18-20, for a demonstration of the stability argument.) 

Whereas the original differential equations of motion are identically 
conservative of energy, the finite difference approximations are not. 

This is an important difference from the situation in Ibrt I, where the 
addition of an internal variable resulted in complete conservation* It 
can te show that In Part II the diserepanoy par cycle la proportional to 
(5t) , so that the cumulative discrepancy over a prescribed elapsed time 
is proportional to 6t and can therefore be made as small as desired. 
Experimentation has shown further that when the overall relative dis¬ 
crepancy in energy has been reduced to a certain smallness, then any addi¬ 
tional decrease in the value of 6t will make a negligible change in any 
of the calculated results, indicating that the solution is as accurate 
as desired. In all calculations used in this study, the value of 6t was 
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indeed sufficiently small that the solutions, in effect, axe those of the 
differential equations. 

The details of the calculations are more particularly defined by spe¬ 
cifying the units of distance and time. (The units of mass do not enter, 
all quantities such as energy and force being measured per unit mass.) 

In all cases, the unit of distance vas chosen to be the initial particle 
spacing, 8 s 1.0. Thus, with N particles, including those at the walls, 
the region of allowed motion was of length (H - 1 )s, and contained (N - 2) 
moving particles. The units of time are completely specified by the 
numerical value of c, which has the dimensions of velocity. For all 
problems, ^ = 1000. The only variables from problem to problem were N 
^hich was either 2? or 52) and U (which was given various positive values 

in the interval 5 < U £ 50) • 

A typical example of the results ohtaincd is given an Fig. II-1 , m 
which N = 27, U = 25. The lines show the positions of all 25 
interior particles as functions of time. It shouM he emphasised that 
the jaggedness of the lines in some parts has resulted from straight-line 
connection of datum points obtained from the computer at sampUng times 
which were many cycles apart. Data plotted from every time cycle would 
Show completely smooth particle trajectories, but were not obtained he- 
cause of the extra computer time required. It should also be noted that 
the periodicity of asmntude oscillation for the particles in the right- 
hand compression region also is an lUnsion resulting from the sUS>t 
difference in phase between the true oscillations and the data 





sampling period. (Samples taken in perfect phase vith the txnie oscilla¬ 
tions wuld have ohsc\ired their actual amplitude, which is a quantity of 
primary interest for the study.) A detailed interpretation of the many 
features of these and similar results forms the main body of this part 
of the report. 

Figure II-2 shows, for the same calculation illustrated in Fig. II-l, 
the histories of total kinetic and potential energy. Conservation of 
total energy can be seen to be almost perfect. The almost exact inter¬ 
change of kinetic and potential energies at later times is a feature to 
be noted; with dissipative forces, the kinetic energy would, instead, be 
e^^cted to be less. 
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Fig. II-2. Kinetic and potential energies as functions of time 
for the calculation illustrated in Fig. II-1. 




CBAHER 2 


MICROSCOPIC UnTERPRETATION 

Our principal purpose in performing the study was to elucidate the 
fluid-dynamics representation of particle-and-force models; nevertheless 
it is helpful to discuss first the results related to molecular dynamics 
inbei^retations. Several of the results will then be useful in the fluid- 

model discussion in Chapter 5 following. 

There are, of course, several respects in which the one-dimensional 

particle array differs from a true molecular assonblage. The primary 
difference is the constraint to one-dimensionality itself; another sig¬ 
nificant difference is the enormous discrepancy between the number of mole- 
cvaes in a nacroaeaplc naterlal swple (~1o“) and the number of partlclee 
in a (~'0^). The first of these differences preelnded the 

possibility of representing a physically realistic problem in these pre¬ 
liminary calculations; the second introduces a considerably closer study 
of fluctiiations than is usually encountered in the kinetic theory of 
gases. The restriction to one-dimensionality nevertheless leaves a number 
of features worth discussing, but any further studies intended to match 
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true moleoJ^ systems must have the restriction removed. Concerning 
the restriction to te^ particles, there already exist extensive data 
£rom some studies by Alder and Wainwri^t [II-4] vbich show that in at 
least one class of problem it is possible to obtain results of consider¬ 
able accuracy vith several hundred particles at most. Their investiga¬ 
tions, however, appear to have been restricted to studies of equilibrium 
properties (or the approach to equilibrium) for particles interacting 
through short-range forces only, and for^ch wall effects have been mini^ 
mized. Our interest, however, is in extending this type of calculation 
to include somewhat lx)nger force ranges in studies of processes far from 
equilibrium. In addition, the effects of walls in studies of shear 
stress, heat conduction, shock generation, etc., are of primary concern, 
and it has been one of our purposes in this and other studies in prog¬ 
ress to see which phenomena couM still be investigated with the limited 
numbers of paartlcles mamgeable by present computers. 

Usually the primary objective in theoretical studies of molecular 
dynamics is to determine the manner in which microscopic phenomena mani¬ 
fest themselves in macroscopic properties. What, for exanrple, are the 
relationships between velocity fluctuations and interparticle potential 
on the one hand to temperature, pressure, and entropy on the other? 

What are the microscopic processes that contribute to viscosity, heat 
conduction, and shock structure? There are, of course, many satisfactory 
answers to these questions for the relatively simple circumstances in 
which inteimolecular forces are weak and departures from equilibrium are 
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small. Beyond these conditions lie the regions where perturbation 
theories no longer apply and the mathematical difficulties seemingly 
become insurmountable. These, then, are the circumstances for which de¬ 
tailed computer solutions of molecular dynamics problems could serve two 
useful purposes. First, they could e:xplain or predict directly the phys¬ 
ical properties of materials subjected to extreme environments. Second, 
they could supply the theoretician with exact and detailed data for com¬ 
parisons with his analysis, for circumstances which are difficult to 
measure e^erimentally • 

The present numerical studies certainly satisfy the conditions of 
being far fran the circumstances in which simple perturbation methods can 
be used for theoretical investigation. For one thing, the intermolecular 
forces are strong, as is demonstrated by the fact that at times the poten¬ 
tial energy of the system significantly exceeds the kinetic energy (see 
Fig. H-2). ]^r another, the processes are strongly time-dependent, in 

volving, in fact, one of the extremes of time-dependent, nonequilibrium 
processes - a strong shock. Thus, we must, for the present, be re¬ 
stricted only to discussions of interpretation of the results, entering 
into analysis for predictions only in heuristic fashion. 

Consider, for example, the shock transition occurring as the par¬ 
ticles "pile ^P" against the right-hand wall (see Fig. H-l). The calcu¬ 
lations ^ow several iniportant features, all consistent with the true pro 
parties of a fluid in such circumstances. Especially evident is the well- 
defined shock itself, which, in its leftward propagation, maintains a 
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sharp transition between two regions whose properties are separately quite 
unifam. The region behind the shock (to the right of it) has properties 
Of particular significance. First, there is a mean compression among the 
particles; second, the mean velocity of the particles has dropped to zero. 
Thus the kinetic energy of mean motion has vanished, and in its place is 
the potential energy acquired by the compression. Most significant of 
all, however, are the residual fluctuations in the compressed region, and 
it is the meaning of these which is worth investigating in detail. 

First, however, note the contrast with the e:^sion region occurring 
at the left side where the particles are moving away from the wall. Again 
there are two regions, each with considerable uniformity, but the transi¬ 
tion is not at all sharp, and, in addition, the stagnant region has no 
fluctuations. 

The detailed e:gjlanation of these phenomena in terms of fluid-dynamics 
interpretations is reserved for Chapter 3i for the present discussion we 
are guided by the fact that the entropy of a fluid increases across a 
shock but remains constant across a rarefaction, together with the infor¬ 
mation that the magnitude of entropy change across a shock can be pre¬ 
dicted from the laws of thermodynamics and the conservation laws of mass, 
nmentum, and energy only, and is ind^endent of any detailed flriid pro¬ 
perties other than those which describe the states on the two sides of 
the shock. It therefore becomes reasonable to inquire whether a detailed 
association can be made between fluctuations and entropy, and there are 
several ways in which such an association can be investigated. Consider 
first the following intuitive argument based upon sinple physical 
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considerations. It has "been shown that the differential equations 
governing the particle motions are completely conservative of energy. If, 
then, there were no fluctuations behind the shock, all of the kinetic 
energy would have been converted to potential, and a subsequent slow ex¬ 
pansion would allow complete recovery of the initial conditions, in viola¬ 
tion of the iireversibility of entropy production throu^ a shock. Thus, 
the shock transition, which results in vanishing kinetic aiergy of mean 
motion, must result in a residual kinetic energy of "random" motion with 
consequent lowering of the recoverable potential energy. Then a subse- 
quenc es^anslon (such as is shown in Fig. II-1 when the left-end rarefac¬ 
tion interacts with the shock) will retxam the system to a state which 
differs fran that at initial time. This is well shown by the retention 
of the fluctuations in Fig. II-l at t « 0.7 when the re-expansion has re¬ 
turned the system nearly to its initial compression; the kinetic energy 
retained in those fluctmtions is that which cannot be recovered, in 
accord with the second law of thermodynamics. 

Yet to be explained, however, is the question of how the paiticles 
"know" just how much fluctuation to create when the shock goes by. As is 
stated before, this cannot depend upon the detailed nat\ire of the inter- 
molecular forces except insofar as they contribute to the pressure states 
on the two sides of the shock, i.e., to the time averages of the inter- 
paiticle forces at only two mean states. This implies that any type of 
detailed variations in intermolecular force which does not change the two 
crucial average forces must result in the same magnlt\jde of shocked-region 
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of the shock transition itself 


fluctuations. Only the detailed structure 
can vary. 

Actually, it is possible to proceed somewhat further in the e^^jlana- 
tion of these phenomena through a consideration of several simple quanti¬ 
tative arguments. We have seen that the equation of motion for particle 
can be written in the form 


2 

d X. 


'IT ' 

dt ^ 


-F(: 



In rtilch F Is a force function ehloh is repulslye nhen posltlre. As a 
model of the behavior of a particle during shock passage, consider the 
motion of when the motions of its neighbors are prescribed to follow 
the me^s e3q)ected for them in the transition, i.e., 

“(jtl)s+ Ut, for t < 

= (J ± l)s + for t > tj^^. 


in which t^^^ is the time at which the shock intersects particle ± l, 
We now define e(t) to represent the deviation of the position of particle 
from the mean position between its two neighbors: 


€(t) s ^x. 



so that e has initial conditions (at t = tj^^) of the form € = 0, dc/dt 
^ Then, except at the two singular times, 
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'"‘J+I ■ ’‘•I-’’ ■ ’ 

dt 

which, for sufficiently small values of €, can he approximated by 

^ = 26 DP ~ s -C^e. 

dt 

Here D symbolizes derivative with respect to argument. The condition 
that is a slowly varying function of time can reasonably be assumed 
in consistency with the smallness of e. [Both require, for validity, 
that U(tj_, - « sj i.e., that the shock be weak.] With slow 

variatiors of A the equation can be solved with sufficiently accurate 
approximation to give 

6(t) = ^ Sin C(t - 


or, for t > tj_^> 


£(t) . ^ .11 ?,*. 

Where we have, without loss of generality, chosen t^^^ = 
^ is given by 


(II-6) 

The quantity 




2 

1 


s -2 



(II-7) 


where 1) is the shock cowresston, avalleble theoretically from the coa 
slderations of Chapter 3. iW-Ued to the force function at hand. 
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F(0 s 


2 


c 


c 


y 


we get 


and 



y 



(II-8) 


Thus the results of this simple analysis are In agreement ulth the 
fact that the ampUtuie of oscillations behind the shock should be Inde¬ 
pendent Of the details of the force function. Furthermore, uhen applied 
to the numerical calculations, the result Is In good agreement. The ob- 
eerved amplitude for the calculation sho™ In Pig. n-,, ^3 

about 0.15. For that problem, U = 25, t = 2.0, c = 31.62, and the ampU- 
tude calculated from Eq. (H-S) is 0.14. 

Prom Eq. (li-6) we see that the fluctuation velocity after shock 
passage is given by 


u 


J 


dt 


I cos ^t. 


(II-9) 


so that the magnitude of the velocity amplitude is indepenaent of the 
force function altogether. Thus the time average of l/2(Uj)^, the 
specific kinetic energy of the fluctuations behind the shock, is 



(II-IO) 


which is one-ei^th of the 


specific kinetic energy in ftront of the shock. 
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Again this result can he compared with those of the calculation shown in 
Fig. II-1, although the comparison cannot he exact, since there is no 
time at which only shocked particles contribute kinetic energy. The nearest 
time mi^ he t = 0.4, at which time nineteen particles have been shocked 
and the others are nearly at rest. Equation CeI- 8) would then predict a 
'total specific kinetic energy 

^ = 742. 

Reference to Pig. II-2, howerer, rtiops the kinetic energy to he somewhat 
than this, and the reason for the discrepancy can he seen in 
Pig. n-1, which shows that particles #6 to 11 have already acquired 
extra hackwards motion from the rarefaction, thereby increasing their 

kinetic energies. 

E:q)loiting even further the results of our simple model, we may de¬ 
rive an expression for the entropy of the assemblage by use of the H 
theorem of Boltzmann, which states that the entropy is given by 

nS = -k // Y In du dx, 

vhere S is the specific entropy of n particles in a region of length L, 

Y is the total one-particle distribution function for the systwn, and 
the integrations extend over all possible positions and velocities. 
Neglecting deviations near walls, I is essentially independent of posi 
tion in a region of uniformity; it represents the probability density 
for finding anjr particle with specified values of position and velocity. 
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and has the normalization property 


// S' du dx = n. 

let u» be the amplitude of velocity fluctuations. Then we may write 

* = (5-)' 

in nhlch (J Is independent of n end L, and depends upon u> only through 

the ratio u/u . Insertion of this into the Boltzmann formula to 

the result 

S s= k ^In ~ + constant^ , 

This, now, can be put into a more useful form throu^ combination with 
Eqs. (II-8) and (II-9), leading finally to 

S = + k In (—) , , 

o Vs /' (n-v 

in which e* is the amplitude of oscillation and is a constant, We 
have therefore verified the correspondence between entropy and oscilla¬ 
tions, showing, in fact, that the former is a function only of the 
amplitude of the latter, as expected. 

The interpretation up to this point is "microscopic," because it is 
based upon the dynamics of the individual particles. Having shown on 
this basis that the microscopic properties of the fluctuations are re¬ 
lated to a gross property (entropy), it will now be possible to proceed 
with the macroscopic interpretation in CSiapter 3. In this 
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interpretation, fluctuation effects are treated as if a dissipative 
force had been added to convert fluctuation energy into the heat of 
macroscopic elements, as would he the case in any realistic general use 
as a fluid-dynamics model. 

There are several additional features that could he studied from 
a molecular-dynamics viewpoint. Some of these have been presented by 
Blackman [II-5] and Butler [II-6], and an extensive report by Gentry, 
Harlw, .and Martin [II-7] vlH soon be published. 
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CHAPTER 3 


MACROSCOPIC INTERHIEEATION 

The investigation presented in Chapter 2 reveals that a nvunber of 
the coorputed particle dynamics features are related to true molecular 
assemblages, and suggests that the statistical pix)pertles of nonequlUbrlum 
processes could be studied vith considerable accuracy through molecular 
trajectory ealcxilations. Such calculations can also be given a different 
interpretation, in idiich each particle r^resents a macroscc^lc element 
of fluid containing an extremely large number of molecules. In reality, 
each element should have an internal degree of freedom representing its 
heat energy. For adiabatic motions, this could be identified with the 
interparticle potential energy, vhile for more general motions, we may 

identify the kinetic energy of fluctuations as the additional dissipative 
heat. 

In the formulation of a practical method, an internal variable is 
introduced for each particle (see Bart I). Then the potential energy is 
fed into this internal variable throu* vork done by the nondlsslpatlve 
force function, vMle a fictitious dissipative force function is intro¬ 
duced to convert the fluctuation energy into addltlorml heat. 
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This aisslpative force snist he veJocity-derpenaent, eo «>et Its effects 
axe irreeerslble end result in removal of the flnctuetlons In vetocity 
(see Bef. n-3). Even vlthout the introduction of this more convenient 
representation of heat energy, it is still possible to give a full, sat¬ 
isfactory interpretation of the present results in terms of fluid-dynamics 
concepts, thereby establishing the basis for development of a practical 

compu'tlBg method. 

One of the basic properties characterising a fluid is its possession 
of an equation of state. This means that the dynamics of a fluid depend 
upon only a few average properties, rather than <®on the detailed motions 
Of all its molecules. Specifically, we mean by equation of state a 
unique relationship among density, heat enersr, and pressure. If such 
can be found for the numerical calculations described in Chapter 2, and 
if its use with fluid-dynamics theory gives results in agreement with 
the macroscopic dynamics observed from the computer, then our demonstra¬ 
tion of a fluid-dynamics interpretation will be complete. This will be 
accodplished in this chapter. 

First, we require an interpretation for density, which has been im¬ 
plied through the introduction of the compression, q, in Chapter 2. 
Worhlng again with specific quantities, the density, p, is defined as 
the reciprocal particle spacing. (The meaning could be made precise 
through calculation of moments of a IdouviUe equation. For our present 
heuristic purposes, the looser definition will be sufficient.) 

Second, we need an erpresslon for heat (internal) energy. 
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TMb comes from Eq. (11.3), In »Mch ire substitute u = u + 8u for the 

velocity, 8u being the fluctuation. Then the local energy per unit mass 
is, on the average. 


E = c ln~ + ~u^+l (6u)^ 

Hov the term u®/2 represents the local klneUc energy of mean motion; all 
the rest must be inteipretea as heat ena^gy. Combining the results of 
Eqs, (II-8) and tn-io;, we find that the eiqnresslon for fluctuation 
kinetic energy can be written ch^/2, m which A - e-,/2. (Thus A Is the 
ratio Of fluctuation amplitude to ^ocal particle spacing.) The required 
expression for specific Internal energy, 1, is thus 


1 . c= in -e. + i c^A® 

Po 2 = * • 


( 11 - 12 ) 


MnalOy, an e^^esslon for pressure is required. This can be derived 
atrectly from the vlrlal theorem In the case of our particular Interparticle 
force function. The details are given by Butler (11-6), tto shows that 


P 


C% (1 + A^). 




[A strai^tforward derivation using the technique of Eef. II-3 yields the 
first part of Eq. (II-13), in confirmtion of the validity of that tech¬ 
nique. Derivation of the second part could likewise be accomplished if a 
proper dissipative force function were chosen; indeed the analysis in that 
case actually would serve as a guide for the choice.] 

These results now can be combined to give the required equation of 

state 


-158- 




p = c^p + 2pl - 2c^p In “ • 

(Nate the factor 2 in the second term, confirmirg agreement vith the sta¬ 
tistical-mechanical equation for a polytropic gas with one degree of free¬ 
dom.) The applicability of this result can he tested in a number of ways, 
we first investigate whether or not the thermodynamic derivation of an ex¬ 
pression for entropy agrees with the statistical result of Chapter 2. 

The starting point is the equation 


p as 

2 S ' 

0 


in which S is a function of p and I. With Eq. this becomes a 

partial differential equation whose solution is 


S •= function of 


• 


The function here Is arbitrary, but a ebolee can be Bade to give complete 
agreement with the statistical result, Eq. (Il-ll). This constitutes 
the first confirmation of the fluid-dynamics interpretation. 

For further investigation we need an expression for the sound speed 

in the fluid, which in general is given by 


W s* 



^const 


A short calculation then gives for this case 



( 11 - 15 ) 
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Consider now the rarefaction at the left in the calculation shown in 
Pig. II- 1 . The envelope of the rarefaction is formed of two sound signals, 
each moving relative to the fluid with ^eed w = c. These calculated 
bounds are shown in Fig. II-1 by the straight lines passing obliquely 
through X = 0 , and it can be seen that agreement with the actual bounds 
is excellent. In addition, the rieJrt-hand bound can be followed alter 
interaction with the shock. According to the theory of Chapter 2 , 

A = U/(2c?y/2), so that Eq. (II-I 5 ) gives w « 35.1 behind the shock. 

Figure II-l shows the continuation with this sound i^ed beyond the 
shock, and again excellent agreement is apparent. 

One additional comparison can be made for the rarefaction region. 
Application of the theory of characteristics to the solution of the rare¬ 
faction problem shows that 

u - c In ~ 

^o 

is a constant through the rarefaction. For the exasrple shown in Fig. H-i, 
this gives the resxilt for density in the jrarefaction region: 

“ * "o (- = <>•'*53 p„. 

Thus the mean distance between particles after rarefaction passage Should 
be 2.21, again in excellent agreement with the observed separation. 

The dynamics of the shock itself, in particular the variation of 
shock propagation velocity with initial material velocity, is the first 
comparison that gives any indication of significant discr^ancy between 



those of the simple heuiristic model used for 


the computer results and 
analysis* 

The shock-dynamics relationships are derived from four basic state¬ 
ments: the three Rankine-Hugoniot conservation principles and the equa¬ 
tion of state. The first three are rigorously correct, so that any dis¬ 
crepancy in the final comparison must bear dijrectly on the adequacy of 
the equation of state. This raises two important questions, which strike 
deeply into the validity of particle-and-force models for fluid dynamics: 
Are the particle dynamics completely capable of representation by an 
equation of state, insofar as macroscopic studies are concerned? Or is 
it possible to put the burden of discrepancy completely onto the approxi¬ 
mations imrolved in deriving the particular equation of state for these 
studies? These questions cannot at present be answered completely, 
although it seems now that the answer to the first is: The dynamics 
are capable of such a representation; the discrepancy is small in 

many circumstances of interest. But such an answer prompts several new 
questions, and in order to answer them, comparisons are still being made 
between particle-dynamics calculations and a variety of fluid dynamics 
situations for which solutions are known. Incidentally, it is easily 
^own that these further investigations, to ^ch a heat variable and 
dissipative force are added, must be carried out in at least two space 
dimensions, where enforced ordering among particles no longer occurs, 
and interpenetrations must be anticipated. Dissipative addition in one 
dimension converts the particle dynamics calculations to familiar 





Iagra n g j .a n counputiog methods which are known to Work prc^erly in a wide 
variety of applications. 

We shall consider two e3q)ressions for the shock speed as a function 
of initial material speed. The first is derived using the conservation 
laws and the full equation of state. The other is a fictitious foxm^ for 
comparison, in \)hich there are forced to be no fluctuations after shock 

passage. For this we use the eqiiations of state in which A = Oj that is, 

2 

p - c p. We shall see that the computer results lie between these two, 
suggesting that the fluctuations predicted by the heuristic analysis are 
somewhat too strong. 

The actual derivations^ which follow standard technique^ need not 
be presented here. For the full equation of state, the algebraic solu¬ 
tion must be performed numerically; for the fictitious case, a simple 
closed solution can be found. The results are presented in Fig. II-3, 
together with the results from saumerous computer calculations in •vdiich 
the initial particle speeds were varied. In all cases, shock velocity 
is relative to material ahead, and thus approaches c = 31.62 as the 
material velocity becomes small. The results themselves give no hltxb as 
to whether the burden of discrepancy should be placed on the capability 
of interpretation or upon the heuristic model used for that purpose. 

It is therefore encouraging that the shocks shown in Part I are in good 
agreement with their eixpected behaviors. 

Other features in Pig. n-i have been calculated, but the results 
add no new extensions to the conclusions beyond those which be 
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SHOCK VELOCITY 



Pig. 11 - 3 • 


Shock ^eed as function of material speed. Cooiputer 
results are shown "by x; solid lines are theoretical. 
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obseirved visually. Thus, for example, it can be seen that the second 
shock in the upper left-hand comer results in a further strengthening 
of the oscillations, as e^qpected, but quantitative prediction yields 
only that same qualitative result. 

Thus we come to a restatement of the principal conclusions from 
Part II of this r^ort; 

1. The Pasrfca-Ulam calculations showed fluid-dynamics properties, 
even though they lacked the two desirable features of correspondence and 
dissipation, because of the close relationship between fluctuations and 

entropy and because of the capability of equation-of-state representation 
even with fluctviations. 

2. The improvements to be expected from the addition of dissipative 
forces and from the proper use of force functions can result in a useful 
fluid-dynamics method, even though the correspondence cannot be “instan¬ 
taneous" but must rely on statistical properties. The results of Part I 
further strengthen this conclusion. 

3. The statistical properties of true molecular assemblages are 
likely to be amenable to computer studies with few particles, even for 
strongly nonequilibrium processes involving Important wall interactions. 
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